#delimit;
clear;
set more off;
version 11;
  quietly log;
  local logon = r(status);
 if "`logon'" == "on" {; log close; };
log using NiemanGibler_JOP_demdiff_appendix.log, text replace;

/* Mark Nieman and Doug Gibler */
/* Appendix for Taking Democratic Differences Seriously */
/* Input: NiemanGibler_demdiff_ally.dta, GiblerTir2014, Threat_MI.dta, 
			Alliance_Formation, NiemanGibler_demdiff_trade */
/* Output: alliance_sel.dta, Figure A1, Figure A2., Figure A3, Figure A4, 
			Table A1-A13 */

/* Preamble */
clear;
estimates clear;
set memory 4g;
program drop _all;
clear matrix;
tempfile temp_file alliance_sel;
set seed 17;

	/* Program for split-sample logit */
		#delimit;
		program define spl_lf, rclass;		
			args lnf beta gamma ;
			tempvar rel violate;
				quietly gen double `rel' = 1/(1+exp(-`gamma'));
				quietly gen double `violate' = 1/(1+exp(-`beta'));
			quietly replace `lnf' = ln((1-`rel')+(`rel'*(1-`violate'))) if $ML_y1==0;	
			quietly replace `lnf' = ln((`rel')*(`violate')) if $ML_y1==1;
		end;
		
		/* program for split-sample probit */
		program define sp_prob_lf, rclass;		
			args lnf beta gamma ;
			tempvar rel violate;
				quietly gen double `rel' = normal(`gamma');
				quietly gen double `violate' = normal(`beta');
			quietly replace `lnf' = ln((1-`rel')+(`rel'*(1-`violate'))) if $ML_y1==0;	
			quietly replace `lnf' = ln((`rel')*(`violate')) if $ML_y1==1;
		end;		
		
	/* Calculate Rubin formula for uncertainty */
		capture program drop rubin;
		program define rubin;
			local names: colnames(Q1);
			/* MI param estimates */
			matrix Q = (Q1 + Q2 + Q3 + Q4 + Q5)/5;
			/* MI covariances     */
			matrix W = (W1 + W2 + W3 + W4 + W5)/5;		
			matrix ll = (ll1 + ll2 + ll3 + ll4 + ll5)/5;
			local  k = colsof(Q);
			forvalues i=1/5 {;
			  matrix QQ=Q`i'-Q;
			  if `i'==1 {;
				matrix B=(QQ)'*QQ;
				};
			  else matrix B = B + (QQ)'*QQ;
			};
			/* Covariance adjustment */
			matrix B=B/(5-1);		
			/* Final covariance   */
			matrix V=W+(1+1/5)*B;				
			ereturn post Q V ;
		end;

/* Construction of the Instrumental Variable */
/* Table A.1: */		
	/* Load Fatal MID and predictors among contigious states */
	use GiblerTir2014, clear;
	/* Replicate Gibler and Tir for time period 1919-2001 */
	logit mzfatalmid samemaster peacetrans violenttrans peace_defpact2 either_cowintra
		threat_militarization1 high_threat_terrMID lndyaddur1 mzpceyrs mzpceyrs2 mzpceyrs3;

		
/* Model Fit */		
/* Table A.2 */
	/* Load alliance data */
	use NiemanGibler_demdiff_ally, clear;
		/* Replicate Leeds et al (Model 1) to identify shared observations*/
		#delimit;
		logit violate wcchangedA capcheither20 regchangeeither chTEAltorBmt alformA demdA 
			asymmetry  nomicoop treaty milinst time timesquare timecube
			, robust cluster(atopidphase);
		gen LMV=1 if e(sample)==1; 

		/* Run Relevancy equation to identify shared observations*/
		logit violate TerrThreat2 num_dems_on_border  num_borders  cold 
			Lenergy_pop Oil_exp  maj_pow kgd_riv , cluster(atopidphase) robust;
		gen GT=1 if e(sample)==1;
		matrix B1 = e(b);

		/* Idenity cases in both samples*/
		gen GL=0;
		replace GL=1 if GT==1 & LMV==1;

		/* Rerun on Leeds et al on same sample (GL) (Model 2)*/
		logit violate wcchangedA capcheither20 regchangeeither chTEAltorBmt alformA demdA 
			asymmetry  nomicoop treaty milinst time timesquare timecube if GL==1
			, robust cluster (atopidphase);
		estimates store LMV2;
		matrix B2 = e(b);
		/* generate individual likelihoods */
			gen double lnfj_lmv = .;
			replace  lnfj_lmv = ln(1/(1+exp((_b[wcchangedA]*wcchangedA + _b[capcheither20]*capcheither20 + _b[regchangeeither]*regchangeeither 
							+ _b[chTEAltorBmt]*chTEAltorBmt + _b[alformA]*alformA + _b[demdA]*demdA + _b[asymmetry]*asymmetry 
							+ _b[nomicoop]*nomicoop + _b[treaty]*treaty + _b[milinst]*milinst + _b[time]*time 
							+ _b[timesquare]*timesquare + _b[timecube]*timecube + _b[_cons])))) if GL==1 & violate==0;
			replace  lnfj_lmv = ln(1/(1+exp(-(_b[wcchangedA]*wcchangedA + _b[capcheither20]*capcheither20 + _b[regchangeeither]*regchangeeither 
							+ _b[chTEAltorBmt]*chTEAltorBmt + _b[alformA]*alformA + _b[demdA]*demdA + _b[asymmetry]*asymmetry 
							+ _b[nomicoop]*nomicoop + _b[treaty]*treaty + _b[milinst]*milinst + _b[time]*time 
							+ _b[timesquare]*timesquare + _b[timecube]*timecube + _b[_cons])))) if GL==1 & violate==1;
		/* Estimate Model 3 with average of Territorial Threat draws */					
		ml model lf spl_lf  
			(B: violate = wcchangedA capcheither20 regchangeeither chTEAltorBmt alformA demdA 
			asymmetry  nomicoop treaty milinst time timesquare timecube)
			(G: violate = Threat_combine num_dems_on_border   num_borders  
			 cold Lenergy_pop Oil_exp maj_pow kgd_riv  ) if GL==1 
			, robust cluster(atopidphase) tech(nr);
		ml init B2 B1 , copy ;
		ml search;
		ml maximize, diff;
		/* Generate individual likelihoods */
			/* split-sample */
			gen double G = 1/(1+exp(-([G]_b[Threat_combine]*Threat_combine +[G]_b[num_dems_on_border]*num_dems_on_border 
							+ [G]_b[num_borders]*num_borders 
							+ [G]_b[cold]*cold + [G]_b[Lenergy_pop]*Lenergy_pop
							+ [G]_b[Oil_exp]*Oil_exp +[G]_b[maj_pow]*maj_pow + [G]_b[kgd_riv]*kgd_riv 
							+ [G]_b[_cons]))) if GL==1;
			gen double B = 1/(1+exp(-([B]_b[wcchangedA]*wcchangedA + [B]_b[capcheither20]*capcheither20 + [B]_b[regchangeeither]*regchangeeither 
							+ [B]_b[chTEAltorBmt]*chTEAltorBmt + [B]_b[alformA]*alformA + [B]_b[demdA]*demdA + [B]_b[asymmetry]*asymmetry 
							+ [B]_b[nomicoop]*nomicoop + [B]_b[treaty]*treaty + [B]_b[milinst]*milinst + [B]_b[time]*time 
							+ [B]_b[timesquare]*timesquare + [B]_b[timecube]*timecube + [B]_b[_cons]))) if GL==1;
			gen double lnfj_full = .;
				replace lnfj_full = ln((1-G)+G*(1-B)) if violate==0 & GL==1;
				replace lnfj_full = ln((G*B)) if violate==1 & GL==1;
	/* Model Comparisons*/
		/* Clarke test (with Schwarz correction)*/
		gen double lnfj_full_adj = lnfj_full*(23/(2*6395))*ln(6395);
		gen double lnfj_lmv_adj = lnfj_lmv*(14/(2*6395))*ln(6395);
		signtest lnfj_full_adj=lnfj_lmv_adj;
		

/* Robustness and Validity Checks */
/* Figure A.1 */
	/* Compare kernel density of alliance violators to sample average - Terr Threat */
	twoway kdensity Threat_combine if violate==1 & GL==1, sort lpattern(solid) lwidth(thick) || kdensity Threat_combine if violate==0 & GL==1, sort scheme(s1mono) lpattern(dash)
					ytitle("Kernel Density") xtitle("Territorial Threat") 
					note("Note: The mean territorial threat for alliance violators is 0.218 with a standard deviation of 0.190,"
						 "while the mean for non-violators is 0.108 with a standard deviation of 0.125. A difference of means"
						 "between the two samples statistically significant at p < 0.001.")
					legend(lab(1  "Violators") lab(2 "Non-violators") symx(*.5));
	graph save figA1, replace asis;
	
/* Table A.3*/	
	/* Note: Relevance is calculated by looking at predicted probabilities using 
	only the relevancy equation ("G" from the split-sample portion of the individual likelihoods) */
	/* Top 15 state-year observations of relevance from sample */
	keep if GL==1;
	collapse (max) Threat_combine G, by(ccode1 year);	
	sort G;
	list ccode1 year G in -15/L;	

/* Table A.4 */
	/* Top 15 state-year observations of threat from sample */
	sort Threat_combine;
	list ccode1 year Threat_combine in -15/L;
	
/* Table A.5 */
	/* Load data set*/	
	use NiemanGibler_demdiff_ally, clear;
	/* Simple */
		/* Replicate Leeds et al to identify shared observations*/
		logit violate  wcchangedA demdA time timesquare timecube
			, robust cluster(atopidphase);
		gen LMV=1 if e(sample)==1; 
		/* Run Relevancy equation to identify shared observations*/
		logit violate TerrThreat1  
			, cluster(atopidphase) robust;
		gen GT=1 if e(sample)==1;
		matrix B1 = e(b);
		/*Idenity cases in both samples*/
		gen GL=0;
		replace GL=1 if GT==1 & LMV==1;
		/* Rerun on Leeds et al on same sample (GL) to get starting values */
		logit violate  wcchangedA demdA time timesquare timecube if GL==1
			, robust cluster (atopidphase);
		matrix B2 = e(b);
		
		/* Simple model */
		forvalues i=1/5 {;
		ml model lf spl_lf  
			(B: violate = wcchangedA demdA time timesquare timecube)
			(G: violate = TerrThreat`i' )if GL==1 
			, robust cluster(atopidphase) tech(nr);
		ml init B2 B1 , copy ;
		ml search;
		ml maximize, diff;	
		matrix Q`i' = e(b);
			matrix W`i' = e(V);
			matrix  ll`i' = e(ll);
			nois _dots `i' 0;
		};
		rubin;
			estout , cells(b(star fmt(3)) se(par fmt(3)) ) starlevels(* .05) style(tex);
			matrix rename ll ll_rub_simple;	
		estimates store simple;
		matrix list ll_rub_simple;		
	
		/*Democratic Development*/
		logit violate TerrThreat1 num_dems_on_border num_borders
			, cluster(atopidphase) robust;
		matrix B1 = e(b);
		forvalues i=1/5 {;
		ml model lf spl_lf  
			(B: violate = wcchangedA demdA time timesquare timecube)
			(G: violate = TerrThreat`i' num_dems_on_border num_borders) if GL==1 
			, robust cluster(atopidphase) tech(nr);
		ml init B2 B1 , copy ;
		ml search;
		ml maximize, diff;
		matrix Q`i' = e(b);
			matrix W`i' = e(V);
			matrix  ll`i' = e(ll);
			nois _dots `i' 0;
		};
		rubin;
			estout , cells(b(star fmt(3)) se(par fmt(3)) ) starlevels(* .05) style(tex);
			matrix rename ll ll_rub_demdev;			
		estimates store demdev;
		matrix list ll_rub_demdev;
		
		/* Full Relevance */
		logit violate TerrThreat1 num_dems_on_border   num_borders  
			 cold Lenergy_pop Oil_exp maj_pow kgd_riv
			, cluster(atopidphase) robust;
		matrix B1 = e(b);
		forvalues i=1/5 {;
		ml model lf spl_lf  
			(B: violate = wcchangedA demdA time timesquare timecube)
			(G: violate = TerrThreat`i' num_dems_on_border   num_borders  
			 cold Lenergy_pop Oil_exp maj_pow kgd_riv) if GL==1 
			, robust cluster(atopidphase) tech(nr);
		ml init B2 B1 , copy ;
		ml search;
		ml maximize, diff;		
		matrix Q`i' = e(b);
			matrix W`i' = e(V);
			matrix  ll`i' = e(ll);
			nois _dots `i' 0;
		};
		rubin;
			estout , cells(b(star fmt(3)) se(par fmt(3)) ) starlevels(* .05) style(tex);
			matrix rename ll ll_rub_fullrel;
		estimates store fullrel;
		matrix list ll_rub_fullrel;		
		
		/* Change in Alliance */
		logit violate wcchangedA demdA capcheither20 regchangeeither chTEAltorBmt alformA time timesquare timecube
			, robust cluster(atopidphase);
		matrix B2 = e(b);
		logit violate TerrThreat1 num_dems_on_border   num_borders  
			 cold Lenergy_pop Oil_exp maj_pow kgd_riv
			, cluster(atopidphase) robust;
		matrix B1 = e(b);
		forvalues i=1/5 {;
		ml model lf spl_lf  
			(B: violate = wcchangedA demdA capcheither20 regchangeeither chTEAltorBmt alformA time timesquare timecube)
			(G: violate = TerrThreat`i' num_dems_on_border   num_borders  
			 cold Lenergy_pop Oil_exp maj_pow kgd_riv) if GL==1 
			, robust cluster(atopidphase) tech(nr);
		ml init B2 B1, copy ;
		ml search;
		ml maximize, diff;
		matrix Q`i' = e(b);
			matrix W`i' = e(V);
			matrix  ll`i' = e(ll);
			nois _dots `i' 0;
		};
		rubin;
			estout , cells(b(star fmt(3)) se(par fmt(3)) ) starlevels(* .05) style(tex);
			matrix rename ll ll_rub_chally;
		estimates store chally;
		matrix list ll_rub_chally;		
	
		/* Alliance Terms */
		logit violate wcchangedA demdA asymmetry  nomicoop treaty milinst time timesquare timecube
			, robust cluster(atopidphase);
		matrix B2 = e(b);
		logit violate TerrThreat1 num_dems_on_border   num_borders  
			 cold Lenergy_pop Oil_exp maj_pow kgd_riv
			, cluster(atopidphase) robust;
		matrix B1 = e(b);
		forvalues i=1/5 {;
		ml model lf spl_lf  
			(B: violate = wcchangedA demdA asymmetry  nomicoop treaty milinst time timesquare timecube)
			(G: violate = TerrThreat`i' num_dems_on_border   num_borders  
			 cold Lenergy_pop Oil_exp maj_pow kgd_riv) if GL==1 
			, robust cluster(atopidphase) tech(nr);
		ml init B2 B1, copy ;
		ml search;
		ml maximize, diff;
		matrix Q`i' = e(b);
			matrix W`i' = e(V);
			matrix  ll`i' = e(ll);
			nois _dots `i' 0;
		};
		rubin;
			estout , cells(b(star fmt(3)) se(par fmt(3)) ) starlevels(* .05) style(tex);
			matrix rename ll ll_rub_ally;
		estimates store ally;	
		matrix list ll_rub_ally;		
		
		/* Full Reliability */
		logit violate wcchangedA demdA capcheither20 regchangeeither chTEAltorBmt alformA asymmetry  nomicoop treaty milinst time timesquare timecube
			, robust cluster(atopidphase);
		matrix B2 = e(b);
		logit violate TerrThreat1 num_dems_on_border num_borders cold 
			, cluster(atopidphase) robust;
		matrix B1 = e(b);
		forvalues i=1/5 {;
		ml model lf spl_lf  
			(B: violate = wcchangedA demdA capcheither20 regchangeeither chTEAltorBmt alformA asymmetry  nomicoop treaty milinst time timesquare timecube)
			(G: violate = TerrThreat`i' num_dems_on_border num_borders cold) if GL==1 
			, robust cluster(atopidphase) tech(nr);
		ml init B2 B1, copy ;
		ml search;
		ml maximize, diff;
		matrix Q`i' = e(b);
			matrix W`i' = e(V);
			matrix  ll`i' = e(ll);
			nois _dots `i' 0;
		};
		rubin;
			estout , cells(b(star fmt(3)) se(par fmt(3)) ) starlevels(* .05) style(tex);
			matrix rename ll ll_rub_reliable;
		estimates store reliable;
		matrix list ll_rub_reliable;		
	
	/* Table A.5 Results */
	estout simple demdev fullrel chally ally reliable, stats(ll N)
	cells(b(fmt(3) star) se(par))
	modelwidth(8)
	starlevels(* 0.05) legend
	varwidth(15)
	style(tex) replace;	
	
/* Table A.6 */
		/* Democracy Directly */
		logit violate wcchangedA capcheither20 regchangeeither chTEAltorBmt alformA demdA asymmetry  nomicoop treaty milinst time timesquare timecube
			, robust cluster(atopidphase);
			matrix B2 = e(b);
		logit violate TerrThreat1 demdA cold Lenergy_pop Oil_exp maj_pow kgd_riv
			, robust cluster(atopidphase);
		matrix B1 = e(b);	
		#delimit;
		forvalues i=1/5 {;
		ml model lf spl_lf  
			(B: violate = wcchangedA capcheither20 regchangeeither chTEAltorBmt alformA demdA 
			asymmetry  nomicoop treaty milinst time timesquare timecube)
			(G: violate = TerrThreat`i' demdA     
			 cold Lenergy_pop Oil_exp maj_pow kgd_riv )if GL==1 
			, robust cluster(atopidphase) tech(bhhh) ;
		ml init B2 B1 , copy ;
		ml search;
		ml maximize, diff;	
		matrix Q`i' = e(b);
			matrix W`i' = e(V);
			matrix  ll`i' = e(ll);
			nois _dots `i' 0;
		};
		rubin;
			estout , cells(b(star fmt(3)) se(par fmt(3)) ) starlevels(* .05) style(tex);
			matrix rename ll ll_rub_dem;	
		estimates store dem;
		matrix list ll_rub_dem;		
		
		/* No Democracy */
		logit violate wcchangedA capcheither20 regchangeeither chTEAltorBmt alformA demdA asymmetry  nomicoop treaty milinst time timesquare timecube
			, robust cluster(atopidphase);
			matrix B2 = e(b);
		logit violate TerrThreat1 num_borders cold Lenergy_pop Oil_exp maj_pow kgd_riv
			, robust cluster(atopidphase);
		matrix B1 = e(b);	
		forvalues i=1/5 {;
		ml model lf spl_lf  
			(B: violate = wcchangedA capcheither20 regchangeeither chTEAltorBmt alformA demdA 
			asymmetry  nomicoop treaty milinst time timesquare timecube)
			(G: violate = TerrThreat`i'  num_borders  
			 cold Lenergy_pop Oil_exp maj_pow kgd_riv )if GL==1 
			, robust cluster(atopidphase) tech(bhhh) ;
		ml init B2 B1 , copy ;
		ml search;
		ml maximize, diff;	
		matrix Q`i' = e(b);
			matrix W`i' = e(V);
			matrix  ll`i' = e(ll);
			nois _dots `i' 0;
		};
		rubin;
			estout , cells(b(star fmt(3)) se(par fmt(3)) ) starlevels(* .05) style(tex);
			matrix rename ll ll_rub_nodem;	
		estimates store nodem;	
		matrix list ll_rub_nodem;		
		
		/* Weak link */
		use Threat_MI, clear;
		gen Threat_combine = (TerrThreat1 + TerrThreat2 + TerrThreat3 + TerrThreat4 
			+ TerrThreat5 + TerrThreat6 + TerrThreat7 + TerrThreat8 + TerrThreat9 + TerrThreat10)/10;
		rename Threat_combine Threat_combineB;
		rename ccode1 ccode2;
		merge 1:m ccode2 year using NiemanGibler_demdiff_ally, gen(_merge_threatB);
		keep if _merge_threatB==3;
		/*weak link measure of threat */
		gen threatweak = min(Threat_combine, Threat_combineB);
		/* Replicate Leeds et al to identify shared observations*/
		logit violate wcchangedA capcheither20 regchangeeither chTEAltorBmt alformA demdA 
			asymmetry  nomicoop treaty milinst time timesquare timecube
			, robust cluster(atopidphase);
		gen LMV=1 if e(sample)==1; 
		matrix lmv_b = e(b);
		/* Run Relevancy equation to identify shared observations*/
		logit violate threatweak num_dems_on_border  num_borders  cold 
			Lenergy_pop Oil_exp  maj_pow kgd_riv , cluster(atopidphase) robust;
		matrix B1_weak = e(b);	
		gen GT=1 if e(sample)==1;
		matrix B1 = e(b);
		/*Idenity cases in both samples*/
		gen GL=0;
		replace GL=1 if GT==1 & LMV==1;
		/* Rerun on Leeds et al on same sample (GL) to get starting values */
		logit violate wcchangedA capcheither20 regchangeeither chTEAltorBmt alformA demdA 
			asymmetry  nomicoop treaty milinst time timesquare timecube if GL==1
			, robust cluster (atopidphase);
		matrix B2 = e(b);		
		/* Estimates with weak link of threat */
		ml model lf spl_lf  
			(B: violate = wcchangedA capcheither20 regchangeeither chTEAltorBmt alformA demdA 
			asymmetry  nomicoop treaty milinst time timesquare timecube)
			(G: violate = threatweak num_dems_on_border   num_borders  
			 cold Lenergy_pop Oil_exp maj_pow kgd_riv  ) if GL==1 
			, robust cluster(atopidphase) tech(nr);
		ml init B2 B1_weak , copy ;
		ml search;
		ml maximize, diff;
		estout , cells(b(star fmt(3)) se(par fmt(3)) ) starlevels(* .05 ) style(tex);		
		estimates store ThreatWeak;	
		
		/* Democratic Allies*/	
		/* Load data set*/
		use NiemanGibler_demdiff_ally, clear;
		gen DemAlliance = 0;
			replace DemAlliance=1 if demdA==1 & demdB==1;
		/* Replicate Leeds et al to identify shared observations*/
		logit violate wcchangedA capcheither20 regchangeeither chTEAltorBmt alformA demdA 
			asymmetry  nomicoop treaty milinst time timesquare timecube
			, robust cluster(atopidphase);
		gen LMV=1 if e(sample)==1; 
		/* Run Relevancy equation to identify shared observations*/
		logit violate Threat_combine DemAlliance num_dems_on_border  num_borders  cold 
			Lenergy_pop Oil_exp  maj_pow kgd_riv , cluster(atopidphase) robust;
		gen GT=1 if e(sample)==1;
		matrix B1 = e(b);
		/*Idenity cases in both samples*/
		gen GL=0;
		replace GL=1 if GT==1 & LMV==1;
		/* Rerun on Leeds et al on same sample (GL) to get starting values */
		logit violate wcchangedA capcheither20 regchangeeither chTEAltorBmt alformA demdA 
			asymmetry  nomicoop treaty milinst time timesquare timecube if GL==1
			, robust cluster (atopidphase);
		matrix B2 = e(b);
		/* Estimates with democratic alliances */
		ml model lf spl_lf  
			(B: violate = wcchangedA capcheither20 regchangeeither chTEAltorBmt alformA demdA 
			asymmetry  nomicoop treaty milinst time timesquare timecube)
			(G: violate = Threat_combine DemAlliance num_dems_on_border   num_borders  
			 cold Lenergy_pop Oil_exp maj_pow kgd_riv  ) if GL==1 
			, robust cluster(atopidphase) tech(nr);
		ml init B2 B1 , copy ;
		ml search;
		ml maximize, diff;
		estout , cells(b(star fmt(3)) se(par fmt(3)) ) starlevels(* .05 ) style(tex);
		estimates store DemAllies;	
			
	/* Table A.6 Results */
	estout dem nodem ThreatWeak DemAllies, stats(ll N)
	cells(b(fmt(3) star) se(par))
	modelwidth(8)
	starlevels(* 0.05) legend
	varwidth(15)
	style(tex) replace;				

/* Table A.7 */
		/* Threat Formation */
		/* Load data */
		use NiemanGibler_demdiff_ally, clear;
		/* Generate threat at alliance formation */
		egen dyad = group(ccode1 ccode2 atopidphase);
		sort dyad year;
		by dyad: gen _t = _n;
		by dyad: gen tf1 = TerrThreat1 if _t==1;
		egen Threat_form1 = max(tf1), by(dyad);
		by dyad: gen tf2 = TerrThreat2 if _t==1;
		egen Threat_form2 = max(tf2), by(dyad);
		by dyad: gen tf3 = TerrThreat3 if _t==1;
		egen Threat_form3 = max(tf3), by(dyad);
		by dyad: gen tf4 = TerrThreat4 if _t==1;
		egen Threat_form4 = max(tf4), by(dyad);
		by dyad: gen tf5 = TerrThreat5 if _t==1;
		egen Threat_form5 = max(tf5), by(dyad);
		/* Replicate Leeds et al to identify shared observations*/
		logit violate wcchangedA capcheither20 regchangeeither chTEAltorBmt alformA demdA 
			asymmetry  nomicoop treaty milinst time timesquare timecube
			, robust cluster(atopidphase);
		gen LMV=1 if e(sample)==1; 
		matrix lmv_b = e(b);
		/* Run Relevancy equation to identify shared observations*/
		logit violate Threat_form2 num_dems_on_border  num_borders  cold 
			Lenergy_pop Oil_exp  maj_pow kgd_riv, cluster(atopidphase) robust ;
		gen GT=1 if e(sample)==1;
		/*Idenity cases in both samples*/
		gen GL=0;
		replace GL=1 if GT==1 & LMV==1;
		/* Rerun on Leeds et al on same sample (GL)*/
		logit violate wcchangedA capcheither20 regchangeeither chTEAltorBmt alformA demdA 
			asymmetry  nomicoop treaty milinst time timesquare timecube if GL==1
			, robust cluster (atopidphase);
		matrix B2 = e(b);
		/*get starting values */
		logit violate Threat_form2 num_dems_on_border  num_borders  cold 
			Lenergy_pop Oil_exp  maj_pow kgd_riv if GL==1, cluster(atopidphase) robust ;
		matrix B1 = e(b);		
		/* Estimate with threat formation */
		forvalues i=1/5 {;
		ml model lf spl_lf  
			(B: violate = wcchangedA capcheither20 regchangeeither chTEAltorBmt alformA demdA 
			asymmetry  nomicoop treaty milinst time timesquare timecube)
			(G: violate = Threat_form`i' num_dems_on_border   num_borders  
			 cold Lenergy_pop Oil_exp maj_pow kgd_riv  ) if GL==1 
			, robust cluster(atopidphase) tech(nr);
		ml init B2 B1 , copy ;
		ml search;
		ml maximize, diff;
		matrix Q`i' = e(b);
			matrix W`i' = e(V);
			matrix  ll`i' = e(ll);
			matrix chi`i' = e(chi2);
			nois _dots `i' 0;
		};
		rubin;
		estout , cells(b(star fmt(3)) se(par fmt(3)) ) starlevels(* .05 ) style(tex);
		estimates store threatform; 
		
		/* Threat formation and time-varying threat measure */
		/* Obtain starting values */
		logit violate Threat_form2 TerrThreat2 num_dems_on_border  num_borders  cold 
			Lenergy_pop Oil_exp  maj_pow kgd_riv if GL==1, cluster(atopidphase) robust ;
		matrix B3 = e(b);
		/* Estimate with threat formation and time-varying measure */
		forvalues i=1/5 {;
		ml model lf spl_lf  
			(B: violate = wcchangedA capcheither20 regchangeeither chTEAltorBmt alformA demdA 
			asymmetry  nomicoop treaty milinst time timesquare timecube)
			(G: violate = Threat_form`i' TerrThreat`i' num_dems_on_border   num_borders  
			 cold Lenergy_pop Oil_exp maj_pow kgd_riv  ) if GL==1 
			, robust cluster(atopidphase) tech(nr);
		ml init B2 B3 , copy ;
		ml search;
		ml maximize, diff;
		matrix Q`i' = e(b);
			matrix W`i' = e(V);
			matrix  ll`i' = e(ll);
			matrix chi`i' = e(chi2);
			nois _dots `i' 0;
		};
		rubin;
		estout , cells(b(star fmt(3)) se(par fmt(3)) ) starlevels(* .05 ) style(tex);
		estimates store threatformtime;
		
		/* Neighbor Average (Note: average neighbor is rescaled by multiplying it by 100) */
		/* Load data */
		use Threat_MI, clear;
		merge 1:m ccode1 year using NiemanGibler_demdiff_ally, gen(_merge_threatavg);
		/* Replicate Leeds et al to identify shared observations*/
		logit violate wcchangedA capcheither20 regchangeeither chTEAltorBmt alformA demdA 
			asymmetry  nomicoop treaty milinst time timesquare timecube
			, robust cluster(atopidphase);
		gen LMV=1 if e(sample)==1; 
		/* Run Relevancy equation to identify shared observations*/
		logit violate TerrThreatAvg1 num_dems_on_border  num_borders  cold 
			Lenergy_pop Oil_exp  maj_pow kgd_riv , cluster(atopidphase) robust;
		gen GT=1 if e(sample)==1;
		matrix B1 = e(b);
		/*Idenity cases in both samples*/
		gen GL=0;
		replace GL=1 if GT==1 & LMV==1;
		/*get starting values */		
		logit violate wcchangedA capcheither20 regchangeeither chTEAltorBmt alformA demdA 
			asymmetry  nomicoop treaty milinst time timesquare timecube if GL==1
			, robust cluster (atopidphase);
		matrix B2 = e(b);		
		/* Estimate with neighbor's average */
		forvalues i=1/5 {;
		ml model lf spl_lf  
			(B: violate = wcchangedA capcheither20 regchangeeither chTEAltorBmt alformA demdA 
			asymmetry  nomicoop treaty milinst time timesquare timecube)
			(G: violate = TerrThreatAvg`i' num_dems_on_border   num_borders  
			 cold Lenergy_pop Oil_exp maj_pow kgd_riv  ) if GL==1 
			, robust cluster(atopidphase) tech(nr);
		ml init B2 B1 , copy ;
		ml search;
		ml maximize, diff;
		matrix Q`i' = e(b);
			matrix W`i' = e(V);
			matrix  ll`i' = e(ll);
			matrix chi`i' = e(chi2);
			nois _dots `i' 0;
		};
		rubin;	
			estout , cells(b(star fmt(3)) se(par fmt(3)) ) starlevels(* .05 ) style(tex);
			matrix rename ll ll_avg;		
			matrix list ll_avg;	
			estimates store NeighborAVG;	
			
		/* 1950-2001 Sample */
		use NiemanGibler_demdiff_ally, clear;
		logit violate TerrThreat2 num_dems_on_border  num_borders  cold 
			Lenergy_pop Oil_exp  maj_pow kgd_riv, cluster(atopidphase) robust;
		matrix B1 = e(b);
		/* Run Relevancy equation to identify shared observations*/
		logit violate wcchangedA capcheither20 regchangeeither chTEAltorBmt alformA demdA 
			asymmetry  nomicoop treaty milinst time timesquare timecube
			, robust cluster (atopidphase);
		matrix B2 = e(b);
		/* Estimate 1950-2001 */
		forvalues i=1/5 {;
		ml model lf spl_lf  
			(B: violate = wcchangedA capcheither20 regchangeeither chTEAltorBmt alformA demdA 
			asymmetry  nomicoop treaty milinst time timesquare timecube)
			(G: violate = TerrThreat`i' num_dems_on_border   num_borders  
			 cold Lenergy_pop Oil_exp maj_pow kgd_riv  ) if  year>1949 
			, robust cluster(atopidphase) tech(nr);
		ml init B2 B1 , copy ;
		ml search;
		ml maximize, diff;
		matrix Q`i' = e(b);
			matrix W`i' = e(V);
			matrix  ll`i' = e(ll);
			matrix chi`i' = e(chi2);
			nois _dots `i' 0;
		};
		rubin;
			estout , cells(b(star fmt(3)) se(par fmt(3)) ) starlevels(* .05 ) style(tex);
			matrix rename ll ll_1950;
			svmat ll_1950;
		estimates store sample_1950; 
		
	/* Table A.7 Results */	
	estout threatform threatformtime NeighborAVG sample_1950, stats(ll N)
	cells(b(fmt(3) star) se(par))
	modelwidth(8)
	starlevels(* 0.05) legend
	varwidth(15)
	style(tex) replace;	
	
/* Table A.8 */	
		/* Load data */
		use Alliance_formation, clear;
		/* Lai and Reiter replication */
		probit CoW3allied CoW3alliedlag pol5 poldif jtrel jtlin jteth mid10ab jtenmid jtmdct10 sqrtdist mpctdum lrnnewd;
			estimates store LaiReiter;
		/* Replicate for 1919-2001 sample */
		probit CoW3allied CoW3alliedlag pol5 poldif jtrel jtlin jteth mid10ab jtenmid jtmdct10 sqrtdist mpctdum lrnnewd if year>1918;
			estimates store LaiReiter1919;
		/* Generate predicted probability of an alliance instrument and inverse Mills ratio */
		predict pr_ally;	
		predict xb_ally;
		gen mill_ally = normalden(xb_ally)/normprob(xb_ally);
		
	/* Table A.8 Results */
	estout LaiReiter LaiReiter1919, stats(ll)
	cells(b(star fmt(3)) se(par fmt(3)) ) 
	starlevels(* .05)  
	style(tex) replace;

/* Table A.9 */		
		/* Set up and merge data with primary data (NiemanGibler_demdiff_ally) */
		keep ccode1 ccode2 year allyforma CoW3allied CoW3alliedlag pol5 poldif jtrel jtlin jteth mid10ab jtenmid jtmdct10 sqrtdist mpctdum lrnnewd pr_ally mill_ally;
		save `temp_file', replace;
			rename ccode1 ccode;
			rename ccode2 ccode1;
			rename ccode ccode2;
		append using `temp_file';
		merge 1:m ccode1 ccode2 year using NiemanGibler_demdiff_ally, gen(_merge_ally_form);
			drop _merge;
		keep if year>1918;
			save alliance_sel, replace;
		/* Replicate Leeds et al (Table 1 Model 1) to identify shared observations*/
		logit violate wcchangedA capcheither20 regchangeeither chTEAltorBmt alformA demdA 
			asymmetry  nomicoop treaty milinst time timesquare timecube 
			, robust cluster(atopidphase);
		gen LMV=1 if e(sample)==1; 
		/* Get initial values for relevance equation without pr_ally (logit) */	
		logit violate Threat_combine num_dems_on_border  num_borders  cold 
			Lenergy_pop Oil_exp  maj_pow kgd_riv, cluster(atopidphase) robust;
		gen GT=1 if e(sample)==1;
		matrix B1 = e(b);
		/*Idenity cases in both samples*/
		gen GL=0;
		replace GL=1 if GT==1 & LMV==1;
		/* Get initial values for relevance equation without mill_ally (probit) */	
		probit violate Threat_combine num_dems_on_border  num_borders  cold 
			Lenergy_pop Oil_exp  maj_pow kgd_riv, cluster(atopidphase) robust;
			matrix B1_probit = e(b);
		/* Get initial values for relevance equation with pr_ally instrument */
		logit violate Threat_combine num_dems_on_border  num_borders  cold 
			Lenergy_pop Oil_exp  maj_pow kgd_riv pr_ally, cluster(atopidphase) robust;
		matrix B1_ia = e(b);
		/* Get initial values for relevance equation with mill_ally instrument */
		probit violate Threat_combine num_dems_on_border  num_borders  cold 
			Lenergy_pop Oil_exp  maj_pow kgd_riv mill_ally, cluster(atopidphase) robust;
		matrix B1_ma = e(b);
		/* Rerun on Leeds et al on same sample (GL) - logit starting values */
		#delimit;
		logit violate wcchangedA capcheither20 regchangeeither chTEAltorBmt alformA demdA 
			asymmetry  nomicoop treaty milinst time timesquare timecube if GL==1
			, robust cluster (atopidphase);
		matrix B2 = e(b);
		/* Rerun on Leeds et al on same sample (GL) - probit starting values */
		probit violate wcchangedA capcheither20 regchangeeither chTEAltorBmt alformA demdA 
			asymmetry  nomicoop treaty milinst time timesquare timecube if GL==1
			, robust cluster (atopidphase);
		matrix B2_probit = e(b);
		
		/* split sample logit with instrument (pr_ally) in 1) relevance, 2) outcome, 3) both */
		ml model lf spl_lf  
			(V: violate = wcchangedA capcheither20 regchangeeither chTEAltorBmt alformA demdA 
			asymmetry  nomicoop treaty milinst time timesquare timecube)
			(R: violate = Threat_combine num_dems_on_border   num_borders  
			 cold Lenergy_pop Oil_exp maj_pow kgd_riv  pr_ally) 
			, robust cluster(atopidphase) tech(nr);
		ml init B2 B1_ia , copy ;
		ml search;
		ml maximize, diff;
		estimates store spl_prallyInstR;	
	
		ml model lf spl_lf  
			(V: violate = wcchangedA capcheither20 regchangeeither chTEAltorBmt alformA demdA 
			asymmetry  nomicoop treaty milinst time timesquare timecube pr_ally)
			(R: violate = Threat_combine num_dems_on_border   num_borders  
			 cold Lenergy_pop Oil_exp maj_pow kgd_riv ) 
			, robust cluster(atopidphase) tech(nr);
		ml init B2 0 B1 , copy ;
		ml search;
		ml maximize, diff;
		estimates store spl_prallyInstV;
				
		ml model lf spl_lf  
			(V: violate = wcchangedA capcheither20 regchangeeither chTEAltorBmt alformA demdA 
			asymmetry  nomicoop treaty milinst time timesquare timecube pr_ally)
			(R: violate = Threat_combine num_dems_on_border   num_borders  
			 cold Lenergy_pop Oil_exp maj_pow kgd_riv  pr_ally) 
			, robust cluster(atopidphase) tech(nr);
		ml init B2 0 B1_ia , copy ;
		ml search;
		ml maximize, diff;
		estimates store spl_prallyInstB;

		/*split sample probit with instrument (mill_ally) in 1) relevance, 2) outcome, 3) both*/
		ml model lf sp_prob_lf  
			(V: violate = wcchangedA capcheither20 regchangeeither chTEAltorBmt alformA demdA 
			asymmetry  nomicoop treaty milinst time timesquare timecube )
			(R: violate = Threat_combine num_dems_on_border   num_borders  
			 cold Lenergy_pop Oil_exp maj_pow kgd_riv mill_ally) if GL==1 
			, robust cluster(atopidphase) tech(nr);
		ml init B2_probit B1_ma, copy ;
		ml search;
		ml maximize, diff;	
		estimates store sp_prob_millallyR;
		
		ml model lf sp_prob_lf  
			(V: violate = wcchangedA capcheither20 regchangeeither chTEAltorBmt alformA demdA 
			asymmetry  nomicoop treaty milinst time timesquare timecube mill_ally)
			(R: violate = Threat_combine num_dems_on_border   num_borders  
			 cold Lenergy_pop Oil_exp maj_pow kgd_riv ) if GL==1 
			, robust cluster(atopidphase) tech(nr);
		ml init B2_probit 0 B1, copy ;
		ml search;
		ml maximize, diff;	
		estimates store sp_prob_millallyV;	
		
		ml model lf sp_prob_lf  
			(V: violate = wcchangedA capcheither20 regchangeeither chTEAltorBmt alformA demdA 
			asymmetry  nomicoop treaty milinst time timesquare timecube mill_ally)
			(R: violate = Threat_combine num_dems_on_border   num_borders  
			 cold Lenergy_pop Oil_exp maj_pow kgd_riv mill_ally) if GL==1 
			, robust cluster(atopidphase) tech(nr);
		ml init B2_probit 0 B1_ma, copy ;
		ml search;
		ml maximize, diff;	
		estimates store sp_prob_millallyB;	
		
	/* Table A.9 Results */
	estout sp*, stats(ll)
	cells(b(star fmt(3)) se(par fmt(3)) ) 
	starlevels(* 0.05) legend    
	style(tex);	
	
/* Table A.10 */
		/* Censored probit of alliance formation and violation outcome (omits relevance) */
		heckprob violate wcchangedA capcheither20 regchangeeither chTEAltorBmt alformA demdA asymmetry  nomicoop treaty milinst time timesquare timecube
		 , sel(CoW3allied = CoW3alliedlag pol5 poldif jtrel jtlin jteth mid10ab jtenmid jtmdct10 sqrtdist mpctdum lrnnewd) cluster(ccode1) robust;
			estimates store heck_orig;

		/* Calculate Pr(R_i) to insert in censored probit outcome equation */
		ml model lf sp_prob_lf  
			(B: violate = wcchangedA capcheither20 regchangeeither chTEAltorBmt alformA demdA 
			asymmetry  nomicoop treaty milinst time timesquare timecube)
			(G: violate = Threat_combine num_dems_on_border   num_borders  
			 cold Lenergy_pop Oil_exp maj_pow kgd_riv  ) 
			, robust cluster(atopidphase) tech(nr);
		ml init B2_probit B1_probit , copy ;
		ml search;
		ml maximize, diff;
		gen double rel = normal([G]_b[Threat_combine]*Threat_combine +[G]_b[num_dems_on_border]*num_dems_on_border 
							+ [G]_b[num_borders]*num_borders 
							+ [G]_b[cold]*cold + [G]_b[Lenergy_pop]*Lenergy_pop
							+ [G]_b[Oil_exp]*Oil_exp +[G]_b[maj_pow]*maj_pow + [G]_b[kgd_riv]*kgd_riv 
							+ [G]_b[_cons]);
		/* Estimate censored probit with Pr(R_i) in outcome equation */
		heckprob violate wcchangedA capcheither20 regchangeeither chTEAltorBmt alformA demdA asymmetry  nomicoop treaty milinst time timesquare timecube rel
		 , sel(CoW3allied = CoW3alliedlag pol5 poldif jtrel jtlin jteth mid10ab jtenmid jtmdct10 sqrtdist mpctdum lrnnewd) cluster(ccode1) robust;
			estimates store heck_rel;
		
		/* Estimate censored probit with all relevance equation variables included the outcome equation */
		heckprob violate wcchangedA capcheither20 regchangeeither chTEAltorBmt alformA demdA asymmetry  nomicoop treaty milinst time timesquare timecube 
		Threat_combine num_dems_on_border num_borders  cold Lenergy_pop Oil_exp maj_pow kgd_riv
		 , sel(CoW3allied = CoW3alliedlag pol5 poldif jtrel jtlin jteth mid10ab jtenmid jtmdct10 sqrtdist mpctdum lrnnewd) cluster(ccode1) robust;
			estimates store heck_all;
		
	/* Table A.10 Results */
	estout heck*,  stats(ll)
	cells(b(star fmt(3)) se(par fmt(3)) ) 
	starlevels(* 0.05) legend 
	style(tex);
	
/* Figure A.2 */
		/* Load data */
		use Threat_MI, clear;
		/* Generate difference in value for t-1, t-2, and t-3 for each draw of territorial threat
		and merge with alliance formation data set */
		sort ccode1 year;
		xtset ccode1 year;
		forvalues i=1/10 {;
			gen ChgThreat_y1_d`i' = D1.TerrThreat`i';
			gen ChgThreat_y2_d`i' = D2.TerrThreat`i';
			gen ChgThreat_y3_d`i' = D3.TerrThreat`i';
		};
		merge 1:m ccode1 year using alliance_sel, gen(_merge_inter);
			drop _merge_inter;
		/* Generage democracy variables */
		gen onedem=1 if demdB==1&demdA==0;
			recode onedem .=1 if demdA==1;
			recode onedem .=0;
		gen jointdem=1 if demdB==1&demdA==1;
			recode jointdem .=0;
		/* Set up alliance data for intervention analysis */
		egen alliance_member=group(atop ccode1);
		sort alliance_member year;
		gen firstatopyear=1 if alliance_member~=alliance_member[_n-1]&atop>2000;
		recode firstatopyear .=0;
		recode firstatopyear 1=. if atop==.;
		/* Generate intervention, post-, and pre-intervention data */
		gen Intervention=0;
			replace Intervention=1 if firstatop==1;
		gen Post1=0;
			replace Post1=1 if firstatop[_n-1]==1 & ccode1==ccode1[_n-1];
		gen Post2=0;
			replace Post2=1 if firstatop[_n-2]==1 & ccode1==ccode1[_n-2];
		gen Post3=0;
			replace Post3=1 if firstatop[_n-3]==1 & ccode1==ccode1[_n-3];
		gen Post4=0;
			replace Post4=1 if firstatop[_n-4]==1 & ccode1==ccode1[_n-4];
		gen Post5=0;
			replace Post5=1 if firstatop[_n-5]==1 & ccode1==ccode1[_n-5];
		egen dyad = group(ccode1 ccode2);
		sort dyad year;
		gen Pre1=0;
			replace Pre1=1 if firstatop[_n+1]==1 & dyad==dyad[_n+1] ;
		gen Pre2=0;
			replace Pre2=1 if firstatop[_n+2]==1 & dyad==dyad[_n+2];
		gen Pre3=0;
			replace Pre3=1 if firstatop[_n+3]==1 & dyad==dyad[_n+3] ;
		gen Pre4=0;
			replace Pre4=1 if firstatop[_n+4]==1  & dyad==dyad[_n+4];
		gen Pre5=0;
			replace Pre5=1 if firstatop[_n+5]==1  & dyad==dyad[_n+5];
		/* Organize data for analysis */
		sort ccode1 ccode2 year	;
		order ccode1 ccode2 year atop firstatopyear Intervention Post* Pre* ChgThreat*;
		/*Analyze Threat t-3 for each value for 5 years pre- and post-intervention */
			/* Apply Rubin for uncertainty in estimates */
			forvalues i=1/10 {;
			regress ChgThreat_y3_d`i' onedem jointdem if Post5==1;
			matrix Q`i' = e(b);
				matrix W`i' = e(V);
				matrix  ll`i' = e(ll);
				matrix chi`i' = e(chi2);
				nois _dots `i' 0;
			};
			rubin;
			estout , cells(b(star fmt(3)) se(par fmt(3)) ) starlevels(* .10 ** .05 *** .01) style(tex);
			/* Obtain data for graph plot */
			gen b3_demA_Post5  = _b[onedem];
			gen se3_lo_demA_Post5 = _b[onedem]-1.96*_se[onedem];
			gen se3_hi_demA_Post5 = _b[onedem]+1.96*_se[onedem];
			gen b3_2dem_Post5  = _b[jointdem];
			gen se3_lo_2dem_Post5 = _b[jointdem]-1.96*_se[jointdem];
			gen se3_hi_2dem_Post5 = _b[jointdem]+1.96*_se[jointdem];
			/* Apply Rubin for uncertainty in estimates */
			forvalues i=1/10 {;
			regress ChgThreat_y3_d`i' onedem jointdem if Post4==1;
			matrix Q`i' = e(b);
				matrix W`i' = e(V);
				matrix  ll`i' = e(ll);
				matrix chi`i' = e(chi2);
				nois _dots `i' 0;
			};
			rubin;
			estout , cells(b(star fmt(3)) se(par fmt(3)) ) starlevels(* .10 ** .05 *** .01) style(tex);
			/* Obtain data for graph plot */
			gen b3_demA_Post4  = _b[onedem];
			gen se3_lo_demA_Post4 = _b[onedem]-1.96*_se[onedem];
			gen se3_hi_demA_Post4 = _b[onedem]+1.96*_se[onedem];
			gen b3_2dem_Post4  = _b[jointdem];
			gen se3_lo_2dem_Post4 = _b[jointdem]-1.96*_se[jointdem];
			gen se3_hi_2dem_Post4 = _b[jointdem]+1.96*_se[jointdem];
			/* Apply Rubin for uncertainty in estimates */
			forvalues i=1/10 {;
			regress ChgThreat_y3_d`i' onedem jointdem if Post3==1;
			matrix Q`i' = e(b);
				matrix W`i' = e(V);
				matrix  ll`i' = e(ll);
				matrix chi`i' = e(chi2);
				nois _dots `i' 0;
			};
			rubin;
			estout , cells(b(star fmt(3)) se(par fmt(3)) ) starlevels(* .10 ** .05 *** .01) style(tex);
			/* Obtain data for graph plot */
			gen b3_demA_Post3  = _b[onedem];
			gen se3_lo_demA_Post3 = _b[onedem]-1.96*_se[onedem];
			gen se3_hi_demA_Post3 = _b[onedem]+1.96*_se[onedem];
			gen b3_2dem_Post3  = _b[jointdem];
			gen se3_lo_2dem_Post3 = _b[jointdem]-1.96*_se[jointdem];
			gen se3_hi_2dem_Post3 = _b[jointdem]+1.96*_se[jointdem];
			/* Apply Rubin for uncertainty in estimates */
			forvalues i=1/10 {;
			regress ChgThreat_y3_d`i' onedem jointdem if Post2==1;
			matrix Q`i' = e(b);
				matrix W`i' = e(V);
				matrix  ll`i' = e(ll);
				matrix chi`i' = e(chi2);
				nois _dots `i' 0;
			};
			rubin;
			estout , cells(b(star fmt(3)) se(par fmt(3)) ) starlevels(* .10 ** .05 *** .01) style(tex);
			/* Obtain data for graph plot */
			gen b3_demA_Post2  = _b[onedem];
			gen se3_lo_demA_Post2 = _b[onedem]-1.96*_se[onedem];
			gen se3_hi_demA_Post2 = _b[onedem]+1.96*_se[onedem];
			gen b3_2dem_Post2  = _b[jointdem];
			gen se3_lo_2dem_Post2 = _b[jointdem]-1.96*_se[jointdem];
			gen se3_hi_2dem_Post2 = _b[jointdem]+1.96*_se[jointdem];
			/* Apply Rubin for uncertainty in estimates */
			forvalues i=1/10 {;
			regress ChgThreat_y3_d`i' onedem jointdem if Post1==1;
			matrix Q`i' = e(b);
				matrix W`i' = e(V);
				matrix  ll`i' = e(ll);
				matrix chi`i' = e(chi2);
				nois _dots `i' 0;
			};
			rubin;
			estout , cells(b(star fmt(3)) se(par fmt(3)) ) starlevels(* .10 ** .05 *** .01) style(tex);
			/* Obtain data for graph plot */
			gen b3_demA_Post1  = _b[onedem];
			gen se3_lo_demA_Post1 = _b[onedem]-1.96*_se[onedem];
			gen se3_hi_demA_Post1 = _b[onedem]+1.96*_se[onedem];
			gen b3_2dem_Post1  = _b[jointdem];
			gen se3_lo_2dem_Post1 = _b[jointdem]-1.96*_se[jointdem];
			gen se3_hi_2dem_Post1 = _b[jointdem]+1.96*_se[jointdem];
			/* Apply Rubin for uncertainty in estimates */
			forvalues i=1/10 {;
			regress ChgThreat_y3_d`i' onedem jointdem if Intervention==1;
			matrix Q`i' = e(b);
				matrix W`i' = e(V);
				matrix  ll`i' = e(ll);
				matrix chi`i' = e(chi2);
				nois _dots `i' 0;
			};
			rubin;
			estout , cells(b(star fmt(3)) se(par fmt(3)) ) starlevels(* .10 ** .05 *** .01) style(tex);
			/* Obtain data for graph plot */
			gen b3_demA_Post0  = _b[onedem];
			gen se3_lo_demA_Post0 = _b[onedem]-1.96*_se[onedem];
			gen se3_hi_demA_Post0 = _b[onedem]+1.96*_se[onedem];
			gen b3_2dem_Post0  = _b[jointdem];
			gen se3_lo_2dem_Post0 = _b[jointdem]-1.96*_se[jointdem];
			gen se3_hi_2dem_Post0 = _b[jointdem]+1.96*_se[jointdem];
			/* Apply Rubin for uncertainty in estimates */
			forvalues i=1/10 {;
			regress ChgThreat_y3_d`i' onedem jointdem if Pre1==1;
			matrix Q`i' = e(b);
				matrix W`i' = e(V);
				matrix  ll`i' = e(ll);
				matrix chi`i' = e(chi2);
				nois _dots `i' 0;
			};
			rubin;
			estout , cells(b(star fmt(3)) se(par fmt(3)) ) starlevels(* .10 ** .05 *** .01) style(tex);
			/* Obtain data for graph plot */
			gen b3_demA_Pre1  = _b[onedem];
			gen se3_lo_demA_Pre1 = _b[onedem]-1.96*_se[onedem];
			gen se3_hi_demA_Pre1 = _b[onedem]+1.96*_se[onedem];
			gen b3_2dem_Pre1  = _b[jointdem];
			gen se3_lo_2dem_Pre1 = _b[jointdem]-1.96*_se[jointdem];
			gen se3_hi_2dem_Pre1 = _b[jointdem]+1.96*_se[jointdem];
			/* Apply Rubin for uncertainty in estimates */
			forvalues i=1/10 {;
			regress ChgThreat_y3_d`i' onedem jointdem if Pre2==1;
			matrix Q`i' = e(b);
				matrix W`i' = e(V);
				matrix  ll`i' = e(ll);
				matrix chi`i' = e(chi2);
				nois _dots `i' 0;
			};
			rubin;
			estout , cells(b(star fmt(3)) se(par fmt(3)) ) starlevels(* .10 ** .05 *** .01) style(tex);
			/* Obtain data for graph plot */
			gen b3_demA_Pre2  = _b[onedem];
			gen se3_lo_demA_Pre2 = _b[onedem]-1.96*_se[onedem];
			gen se3_hi_demA_Pre2 = _b[onedem]+1.96*_se[onedem];
			gen b3_2dem_Pre2  = _b[jointdem];
			gen se3_lo_2dem_Pre2 = _b[jointdem]-1.96*_se[jointdem];
			gen se3_hi_2dem_Pre2 = _b[jointdem]+1.96*_se[jointdem];
			/* Apply Rubin for uncertainty in estimates */
			forvalues i=1/10 {;
			regress ChgThreat_y3_d`i' onedem jointdem if Pre3==1;
			matrix Q`i' = e(b);
				matrix W`i' = e(V);
				matrix  ll`i' = e(ll);
				matrix chi`i' = e(chi2);
				nois _dots `i' 0;
			};
			rubin;
			estout , cells(b(star fmt(3)) se(par fmt(3)) ) starlevels(* .10 ** .05 *** .01) style(tex);
			/* Obtain data for graph plot */
			gen b3_demA_Pre3  = _b[onedem];
			gen se3_lo_demA_Pre3 = _b[onedem]-1.96*_se[onedem];
			gen se3_hi_demA_Pre3 = _b[onedem]+1.96*_se[onedem];
			gen b3_2dem_Pre3  = _b[jointdem];
			gen se3_lo_2dem_Pre3 = _b[jointdem]-1.96*_se[jointdem];
			gen se3_hi_2dem_Pre3 = _b[jointdem]+1.96*_se[jointdem];
			/* Apply Rubin for uncertainty in estimates */
			forvalues i=1/10 {;
			regress ChgThreat_y3_d`i' onedem jointdem if Pre4==1;
			matrix Q`i' = e(b);
				matrix W`i' = e(V);
				matrix  ll`i' = e(ll);
				matrix chi`i' = e(chi2);
				nois _dots `i' 0;
			};
			rubin;
			estout , cells(b(star fmt(3)) se(par fmt(3)) ) starlevels(* .10 ** .05 *** .01) style(tex);
			/* Obtain data for graph plot */
			gen b3_demA_Pre4  = _b[onedem];
			gen se3_lo_demA_Pre4 = _b[onedem]-1.96*_se[onedem];
			gen se3_hi_demA_Pre4 = _b[onedem]+1.96*_se[onedem];
			gen b3_2dem_Pre4  = _b[jointdem];
			gen se3_lo_2dem_Pre4 = _b[jointdem]-1.96*_se[jointdem];
			gen se3_hi_2dem_Pre4 = _b[jointdem]+1.96*_se[jointdem];
			/* Apply Rubin for uncertainty in estimates */
			forvalues i=1/10 {;
			regress ChgThreat_y3_d`i' onedem jointdem if Pre5==1;
			matrix Q`i' = e(b);
				matrix W`i' = e(V);
				matrix  ll`i' = e(ll);
				matrix chi`i' = e(chi2);
				nois _dots `i' 0;
			};
			rubin;
			estout , cells(b(star fmt(3)) se(par fmt(3)) ) starlevels(* .10 ** .05 *** .01) style(tex);
			/* Obtain data for graph plot */
			gen b3_demA_Pre5  = _b[onedem];
			gen se3_lo_demA_Pre5 = _b[onedem]-1.96*_se[onedem];
			gen se3_hi_demA_Pre5 = _b[onedem]+1.96*_se[onedem];
			gen b3_2dem_Pre5  = _b[jointdem];
			gen se3_lo_2dem_Pre5 = _b[jointdem]-1.96*_se[jointdem];
			gen se3_hi_2dem_Pre5 = _b[jointdem]+1.96*_se[jointdem];
			/* Generate placement for graph */
			gen Chg3_pre5=1;
			gen Chg3_pre4=2;
			gen Chg3_pre3=3;
			gen Chg3_pre2=4;
			gen Chg3_pre1=5;
			gen Chg3_post0=6;
			gen Chg3_post1=7;
			gen Chg3_post2=8;
			gen Chg3_post3=9;
			gen Chg3_post4=10;
			gen Chg3_post5=11;


		/*Analyze Threat t-2 for each value for 5 years pre- and post-intervention */
			/* Apply Rubin for uncertainty in estimates */
			forvalues i=1/10 {;
			regress ChgThreat_y2_d`i' onedem jointdem if Post5==1;
			matrix Q`i' = e(b);
				matrix W`i' = e(V);
				matrix  ll`i' = e(ll);
				matrix chi`i' = e(chi2);
				nois _dots `i' 0;
			};
			rubin;
			estout , cells(b(star fmt(3)) se(par fmt(3)) ) starlevels(* .10 ** .05 *** .01) style(tex);
			/* Obtain data for graph plot */
			gen b2_demA_Post5  = _b[onedem];
			gen se2_lo_demA_Post5 = _b[onedem]-1.96*_se[onedem];
			gen se2_hi_demA_Post5 = _b[onedem]+1.96*_se[onedem];
			gen b2_2dem_Post5  = _b[jointdem];
			gen se2_lo_2dem_Post5 = _b[jointdem]-1.96*_se[jointdem];
			gen se2_hi_2dem_Post5 = _b[jointdem]+1.96*_se[jointdem];
			/* Apply Rubin for uncertainty in estimates */
			forvalues i=1/10 {;
			regress ChgThreat_y2_d`i' onedem jointdem if Post4==1;
			matrix Q`i' = e(b);
				matrix W`i' = e(V);
				matrix  ll`i' = e(ll);
				matrix chi`i' = e(chi2);
				nois _dots `i' 0;
			};
			rubin;
			estout , cells(b(star fmt(3)) se(par fmt(3)) ) starlevels(* .10 ** .05 *** .01) style(tex);
			/* Obtain data for graph plot */
			gen b2_demA_Post4  = _b[onedem];
			gen se2_lo_demA_Post4 = _b[onedem]-1.96*_se[onedem];
			gen se2_hi_demA_Post4 = _b[onedem]+1.96*_se[onedem];
			gen b2_2dem_Post4  = _b[jointdem];
			gen se2_lo_2dem_Post4 = _b[jointdem]-1.96*_se[jointdem];
			gen se2_hi_2dem_Post4 = _b[jointdem]+1.96*_se[jointdem];
			/* Apply Rubin for uncertainty in estimates */
			forvalues i=1/10 {;
			regress ChgThreat_y2_d`i' onedem jointdem if Post3==1;
			matrix Q`i' = e(b);
				matrix W`i' = e(V);
				matrix  ll`i' = e(ll);
				matrix chi`i' = e(chi2);
				nois _dots `i' 0;
			};
			rubin;
			estout , cells(b(star fmt(3)) se(par fmt(3)) ) starlevels(* .10 ** .05 *** .01) style(tex);
			/* Obtain data for graph plot */
			gen b2_demA_Post3  = _b[onedem];
			gen se2_lo_demA_Post3 = _b[onedem]-1.96*_se[onedem];
			gen se2_hi_demA_Post3 = _b[onedem]+1.96*_se[onedem];
			gen b2_2dem_Post3  = _b[jointdem];
			gen se2_lo_2dem_Post3 = _b[jointdem]-1.96*_se[jointdem];
			gen se2_hi_2dem_Post3 = _b[jointdem]+1.96*_se[jointdem];
			/* Apply Rubin for uncertainty in estimates */
			forvalues i=1/10 {;
			regress ChgThreat_y2_d`i' onedem jointdem if Post2==1;
			matrix Q`i' = e(b);
				matrix W`i' = e(V);
				matrix  ll`i' = e(ll);
				matrix chi`i' = e(chi2);
				nois _dots `i' 0;
			};
			rubin;
			estout , cells(b(star fmt(3)) se(par fmt(3)) ) starlevels(* .10 ** .05 *** .01) style(tex);
			/* Obtain data for graph plot */
			gen b2_demA_Post2  = _b[onedem];
			gen se2_lo_demA_Post2 = _b[onedem]-1.96*_se[onedem];
			gen se2_hi_demA_Post2 = _b[onedem]+1.96*_se[onedem];
			gen b2_2dem_Post2  = _b[jointdem];
			gen se2_lo_2dem_Post2 = _b[jointdem]-1.96*_se[jointdem];
			gen se2_hi_2dem_Post2 = _b[jointdem]+1.96*_se[jointdem];
			/* Apply Rubin for uncertainty in estimates */
			forvalues i=1/10 {;
			regress ChgThreat_y2_d`i' onedem jointdem if Post1==1;
			matrix Q`i' = e(b);
				matrix W`i' = e(V);
				matrix  ll`i' = e(ll);
				matrix chi`i' = e(chi2);
				nois _dots `i' 0;
			};
			rubin;
			estout , cells(b(star fmt(3)) se(par fmt(3)) ) starlevels(* .10 ** .05 *** .01) style(tex);
			/* Obtain data for graph plot */
			gen b2_demA_Post1  = _b[onedem];
			gen se2_lo_demA_Post1 = _b[onedem]-1.96*_se[onedem];
			gen se2_hi_demA_Post1 = _b[onedem]+1.96*_se[onedem];
			gen b2_2dem_Post1  = _b[jointdem];
			gen se2_lo_2dem_Post1 = _b[jointdem]-1.96*_se[jointdem];
			gen se2_hi_2dem_Post1 = _b[jointdem]+1.96*_se[jointdem];
			/* Apply Rubin for uncertainty in estimates */
			forvalues i=1/10 {;
			regress ChgThreat_y2_d`i' onedem jointdem if Intervention==1;
			matrix Q`i' = e(b);
				matrix W`i' = e(V);
				matrix  ll`i' = e(ll);
				matrix chi`i' = e(chi2);
				nois _dots `i' 0;
			};
			rubin;
			estout , cells(b(star fmt(3)) se(par fmt(3)) ) starlevels(* .10 ** .05 *** .01) style(tex);
			/* Obtain data for graph plot */
			gen b2_demA_Post0  = _b[onedem];
			gen se2_lo_demA_Post0 = _b[onedem]-1.96*_se[onedem];
			gen se2_hi_demA_Post0 = _b[onedem]+1.96*_se[onedem];
			gen b2_2dem_Post0  = _b[jointdem];
			gen se2_lo_2dem_Post0 = _b[jointdem]-1.96*_se[jointdem];
			gen se2_hi_2dem_Post0 = _b[jointdem]+1.96*_se[jointdem];
			/* Apply Rubin for uncertainty in estimates */
			forvalues i=1/10 {;
			regress ChgThreat_y2_d`i' onedem jointdem if Pre1==1;
			matrix Q`i' = e(b);
				matrix W`i' = e(V);
				matrix  ll`i' = e(ll);
				matrix chi`i' = e(chi2);
				nois _dots `i' 0;
			};
			rubin;
			estout , cells(b(star fmt(3)) se(par fmt(3)) ) starlevels(* .10 ** .05 *** .01) style(tex);
			/* Obtain data for graph plot */
			gen b2_demA_Pre1  = _b[onedem];
			gen se2_lo_demA_Pre1 = _b[onedem]-1.96*_se[onedem];
			gen se2_hi_demA_Pre1 = _b[onedem]+1.96*_se[onedem];
			gen b2_2dem_Pre1  = _b[jointdem];
			gen se2_lo_2dem_Pre1 = _b[jointdem]-1.96*_se[jointdem];
			gen se2_hi_2dem_Pre1 = _b[jointdem]+1.96*_se[jointdem];
			/* Apply Rubin for uncertainty in estimates */
			forvalues i=1/10 {;
			regress ChgThreat_y2_d`i' onedem jointdem if Pre2==1;
			matrix Q`i' = e(b);
				matrix W`i' = e(V);
				matrix  ll`i' = e(ll);
				matrix chi`i' = e(chi2);
				nois _dots `i' 0;
			};
			rubin;
			estout , cells(b(star fmt(3)) se(par fmt(3)) ) starlevels(* .10 ** .05 *** .01) style(tex);
			/* Obtain data for graph plot */
			gen b2_demA_Pre2  = _b[onedem];
			gen se2_lo_demA_Pre2 = _b[onedem]-1.96*_se[onedem];
			gen se2_hi_demA_Pre2 = _b[onedem]+1.96*_se[onedem];
			gen b2_2dem_Pre2  = _b[jointdem];
			gen se2_lo_2dem_Pre2 = _b[jointdem]-1.96*_se[jointdem];
			gen se2_hi_2dem_Pre2 = _b[jointdem]+1.96*_se[jointdem];
			/* Apply Rubin for uncertainty in estimates */
			forvalues i=1/10 {;
			regress ChgThreat_y2_d`i' onedem jointdem if Pre3==1;
			matrix Q`i' = e(b);
				matrix W`i' = e(V);
				matrix  ll`i' = e(ll);
				matrix chi`i' = e(chi2);
				nois _dots `i' 0;
			};
			rubin;
			estout , cells(b(star fmt(3)) se(par fmt(3)) ) starlevels(* .10 ** .05 *** .01) style(tex);
			/* Obtain data for graph plot */
			gen b2_demA_Pre3  = _b[onedem];
			gen se2_lo_demA_Pre3 = _b[onedem]-1.96*_se[onedem];
			gen se2_hi_demA_Pre3 = _b[onedem]+1.96*_se[onedem];
			gen b2_2dem_Pre3  = _b[jointdem];
			gen se2_lo_2dem_Pre3 = _b[jointdem]-1.96*_se[jointdem];
			gen se2_hi_2dem_Pre3 = _b[jointdem]+1.96*_se[jointdem];
			/* Apply Rubin for uncertainty in estimates */
			forvalues i=1/10 {;
			regress ChgThreat_y2_d`i' onedem jointdem if Pre4==1;
			matrix Q`i' = e(b);
				matrix W`i' = e(V);
				matrix  ll`i' = e(ll);
				matrix chi`i' = e(chi2);
				nois _dots `i' 0;
			};
			rubin;
			estout , cells(b(star fmt(3)) se(par fmt(3)) ) starlevels(* .10 ** .05 *** .01) style(tex);
			/* Obtain data for graph plot */
			gen b2_demA_Pre4  = _b[onedem];
			gen se2_lo_demA_Pre4 = _b[onedem]-1.96*_se[onedem];
			gen se2_hi_demA_Pre4 = _b[onedem]+1.96*_se[onedem];
			gen b2_2dem_Pre4  = _b[jointdem];
			gen se2_lo_2dem_Pre4 = _b[jointdem]-1.96*_se[jointdem];
			gen se2_hi_2dem_Pre4 = _b[jointdem]+1.96*_se[jointdem];
			/* Apply Rubin for uncertainty in estimates */
			forvalues i=1/10 {;
			regress ChgThreat_y2_d`i' onedem jointdem if Pre5==1;
			matrix Q`i' = e(b);
				matrix W`i' = e(V);
				matrix  ll`i' = e(ll);
				matrix chi`i' = e(chi2);
				nois _dots `i' 0;
			};
			rubin;
			/* Obtain data for graph plot */
			estout , cells(b(star fmt(3)) se(par fmt(3)) ) starlevels(* .10 ** .05 *** .01) style(tex);
			gen b2_demA_Pre5  = _b[onedem];
			gen se2_lo_demA_Pre5 = _b[onedem]-1.96*_se[onedem];
			gen se2_hi_demA_Pre5 = _b[onedem]+1.96*_se[onedem];
			gen b2_2dem_Pre5  = _b[jointdem];
			gen se2_lo_2dem_Pre5 = _b[jointdem]-1.96*_se[jointdem];
			gen se2_hi_2dem_Pre5 = _b[jointdem]+1.96*_se[jointdem];
			/* Generate placement for graph */
			gen Chg2_pre5=1;
			gen Chg2_pre4=2;
			gen Chg2_pre3=3;
			gen Chg2_pre2=4;
			gen Chg2_pre1=5;
			gen Chg2_post0=6;
			gen Chg2_post1=7;
			gen Chg2_post2=8;
			gen Chg2_post3=9;
			gen Chg2_post4=10;
			gen Chg2_post5=11;


		/*Analyze Threat t-1 for each value for 5 years pre- and post-intervention */
			/* Apply Rubin for uncertainty in estimates */
			forvalues i=1/10 {;
			regress ChgThreat_y1_d`i' onedem jointdem if Post5==1;
			matrix Q`i' = e(b);
				matrix W`i' = e(V);
				matrix  ll`i' = e(ll);
				matrix chi`i' = e(chi2);
				nois _dots `i' 0;
			};
			rubin;
			estout , cells(b(star fmt(3)) se(par fmt(3)) ) starlevels(* .10 ** .05 *** .01) style(tex);
			/* Obtain data for graph plot */
			gen b1_demA_Post5  = _b[onedem];
			gen se1_lo_demA_Post5 = _b[onedem]-1.96*_se[onedem];
			gen se1_hi_demA_Post5 = _b[onedem]+1.96*_se[onedem];
			gen b1_2dem_Post5  = _b[jointdem];
			gen se1_lo_2dem_Post5 = _b[jointdem]-1.96*_se[jointdem];
			gen se1_hi_2dem_Post5 = _b[jointdem]+1.96*_se[jointdem];
			/* Apply Rubin for uncertainty in estimates */
			forvalues i=1/10 {;
			regress ChgThreat_y1_d`i' onedem jointdem if Post4==1;
			matrix Q`i' = e(b);
				matrix W`i' = e(V);
				matrix  ll`i' = e(ll);
				matrix chi`i' = e(chi2);
				nois _dots `i' 0;
			};
			rubin;
			estout , cells(b(star fmt(3)) se(par fmt(3)) ) starlevels(* .10 ** .05 *** .01) style(tex);
			/* Obtain data for graph plot */
			gen b1_demA_Post4  = _b[onedem];
			gen se1_lo_demA_Post4 = _b[onedem]-1.96*_se[onedem];
			gen se1_hi_demA_Post4 = _b[onedem]+1.96*_se[onedem];
			gen b1_2dem_Post4  = _b[jointdem];
			gen se1_lo_2dem_Post4 = _b[jointdem]-1.96*_se[jointdem];
			gen se1_hi_2dem_Post4 = _b[jointdem]+1.96*_se[jointdem];
			/* Apply Rubin for uncertainty in estimates */
			forvalues i=1/10 {;
			regress ChgThreat_y1_d`i' onedem jointdem if Post3==1;
			matrix Q`i' = e(b);
				matrix W`i' = e(V);
				matrix  ll`i' = e(ll);
				matrix chi`i' = e(chi2);
				nois _dots `i' 0;
			};
			rubin;
			estout , cells(b(star fmt(3)) se(par fmt(3)) ) starlevels(* .10 ** .05 *** .01) style(tex);
			/* Obtain data for graph plot */
			gen b1_demA_Post3  = _b[onedem];
			gen se1_lo_demA_Post3 = _b[onedem]-1.96*_se[onedem];
			gen se1_hi_demA_Post3 = _b[onedem]+1.96*_se[onedem];
			gen b1_2dem_Post3  = _b[jointdem];
			gen se1_lo_2dem_Post3 = _b[jointdem]-1.96*_se[jointdem];
			gen se1_hi_2dem_Post3 = _b[jointdem]+1.96*_se[jointdem];
			/* Apply Rubin for uncertainty in estimates */
			forvalues i=1/10 {;
			regress ChgThreat_y1_d`i' onedem jointdem if Post2==1;
			matrix Q`i' = e(b);
				matrix W`i' = e(V);
				matrix  ll`i' = e(ll);
				matrix chi`i' = e(chi2);
				nois _dots `i' 0;
			};
			rubin;
			estout , cells(b(star fmt(3)) se(par fmt(3)) ) starlevels(* .10 ** .05 *** .01) style(tex);
			/* Obtain data for graph plot */
			gen b1_demA_Post2  = _b[onedem];
			gen se1_lo_demA_Post2 = _b[onedem]-1.96*_se[onedem];
			gen se1_hi_demA_Post2 = _b[onedem]+1.96*_se[onedem];
			gen b1_2dem_Post2  = _b[jointdem];
			gen se1_lo_2dem_Post2 = _b[jointdem]-1.96*_se[jointdem];
			gen se1_hi_2dem_Post2 = _b[jointdem]+1.96*_se[jointdem];
			/* Apply Rubin for uncertainty in estimates */
			forvalues i=1/10 {;
			regress ChgThreat_y1_d`i' onedem jointdem if Post1==1;
			matrix Q`i' = e(b);
				matrix W`i' = e(V);
				matrix  ll`i' = e(ll);
				matrix chi`i' = e(chi2);
				nois _dots `i' 0;
			};
			rubin;
			estout , cells(b(star fmt(3)) se(par fmt(3)) ) starlevels(* .10 ** .05 *** .01) style(tex);
			/* Obtain data for graph plot */
			gen b1_demA_Post1  = _b[onedem];
			gen se1_lo_demA_Post1 = _b[onedem]-1.96*_se[onedem];
			gen se1_hi_demA_Post1 = _b[onedem]+1.96*_se[onedem];
			gen b1_2dem_Post1  = _b[jointdem];
			gen se1_lo_2dem_Post1 = _b[jointdem]-1.96*_se[jointdem];
			gen se1_hi_2dem_Post1 = _b[jointdem]+1.96*_se[jointdem];
			/* Apply Rubin for uncertainty in estimates */
			forvalues i=1/10 {;
			regress ChgThreat_y1_d`i' onedem jointdem if Intervention==1;
			matrix Q`i' = e(b);
				matrix W`i' = e(V);
				matrix  ll`i' = e(ll);
				matrix chi`i' = e(chi2);
				nois _dots `i' 0;
			};
			rubin;
			estout , cells(b(star fmt(3)) se(par fmt(3)) ) starlevels(* .10 ** .05 *** .01) style(tex);
			/* Obtain data for graph plot */
			gen b1_demA_Post0  = _b[onedem];
			gen se1_lo_demA_Post0 = _b[onedem]-1.96*_se[onedem];
			gen se1_hi_demA_Post0 = _b[onedem]+1.96*_se[onedem];
			gen b1_2dem_Post0  = _b[jointdem];
			gen se1_lo_2dem_Post0 = _b[jointdem]-1.96*_se[jointdem];
			gen se1_hi_2dem_Post0 = _b[jointdem]+1.96*_se[jointdem];
			/* Apply Rubin for uncertainty in estimates */
			forvalues i=1/10 {;
			regress ChgThreat_y1_d`i' onedem jointdem if Pre1==1;
			matrix Q`i' = e(b);
				matrix W`i' = e(V);
				matrix  ll`i' = e(ll);
				matrix chi`i' = e(chi2);
				nois _dots `i' 0;
			};
			rubin;
			estout , cells(b(star fmt(3)) se(par fmt(3)) ) starlevels(* .10 ** .05 *** .01) style(tex);
			/* Obtain data for graph plot */
			gen b1_demA_Pre1  = _b[onedem];
			gen se1_lo_demA_Pre1 = _b[onedem]-1.96*_se[onedem];
			gen se1_hi_demA_Pre1 = _b[onedem]+1.96*_se[onedem];
			gen b1_2dem_Pre1  = _b[jointdem];
			gen se1_lo_2dem_Pre1 = _b[jointdem]-1.96*_se[jointdem];
			gen se1_hi_2dem_Pre1 = _b[jointdem]+1.96*_se[jointdem];
			/* Apply Rubin for uncertainty in estimates */
			forvalues i=1/10 {;
			regress ChgThreat_y1_d`i' onedem jointdem if Pre2==1;
			matrix Q`i' = e(b);
				matrix W`i' = e(V);
				matrix  ll`i' = e(ll);
				matrix chi`i' = e(chi2);
				nois _dots `i' 0;
			};
			rubin;
			estout , cells(b(star fmt(3)) se(par fmt(3)) ) starlevels(* .10 ** .05 *** .01) style(tex);
			/* Obtain data for graph plot */
			gen b1_demA_Pre2  = _b[onedem];
			gen se1_lo_demA_Pre2 = _b[onedem]-1.96*_se[onedem];
			gen se1_hi_demA_Pre2 = _b[onedem]+1.96*_se[onedem];
			gen b1_2dem_Pre2  = _b[jointdem];
			gen se1_lo_2dem_Pre2 = _b[jointdem]-1.96*_se[jointdem];
			gen se1_hi_2dem_Pre2 = _b[jointdem]+1.96*_se[jointdem];
			/* Apply Rubin for uncertainty in estimates */
			forvalues i=1/10 {;
			regress ChgThreat_y1_d`i' onedem jointdem if Pre3==1;
			matrix Q`i' = e(b);
				matrix W`i' = e(V);
				matrix  ll`i' = e(ll);
				matrix chi`i' = e(chi2);
				nois _dots `i' 0;
			};
			rubin;
			estout , cells(b(star fmt(3)) se(par fmt(3)) ) starlevels(* .10 ** .05 *** .01) style(tex);
			/* Obtain data for graph plot */
			gen b1_demA_Pre3  = _b[onedem];
			gen se1_lo_demA_Pre3 = _b[onedem]-1.96*_se[onedem];
			gen se1_hi_demA_Pre3 = _b[onedem]+1.96*_se[onedem];
			gen b1_2dem_Pre3  = _b[jointdem];
			gen se1_lo_2dem_Pre3 = _b[jointdem]-1.96*_se[jointdem];
			gen se1_hi_2dem_Pre3 = _b[jointdem]+1.96*_se[jointdem];
			/* Apply Rubin for uncertainty in estimates */
			forvalues i=1/10 {;
			regress ChgThreat_y1_d`i' onedem jointdem if Pre4==1;
			matrix Q`i' = e(b);
				matrix W`i' = e(V);
				matrix  ll`i' = e(ll);
				matrix chi`i' = e(chi2);
				nois _dots `i' 0;
			};
			rubin;
			estout , cells(b(star fmt(3)) se(par fmt(3)) ) starlevels(* .10 ** .05 *** .01) style(tex);
			/* Obtain data for graph plot */
			gen b1_demA_Pre4  = _b[onedem];
			gen se1_lo_demA_Pre4 = _b[onedem]-1.96*_se[onedem];
			gen se1_hi_demA_Pre4 = _b[onedem]+1.96*_se[onedem];
			gen b1_2dem_Pre4  = _b[jointdem];
			gen se1_lo_2dem_Pre4 = _b[jointdem]-1.96*_se[jointdem];
			gen se1_hi_2dem_Pre4 = _b[jointdem]+1.96*_se[jointdem];
			/* Apply Rubin for uncertainty in estimates */
			forvalues i=1/10 {;
			regress ChgThreat_y1_d`i' onedem jointdem if Pre5==1;
			matrix Q`i' = e(b);
				matrix W`i' = e(V);
				matrix  ll`i' = e(ll);
				matrix chi`i' = e(chi2);
				nois _dots `i' 0;
			};
			rubin;
			estout , cells(b(star fmt(3)) se(par fmt(3)) ) starlevels(* .10 ** .05 *** .01) style(tex);
			/* Obtain data for graph plot */
			gen b1_demA_Pre5  = _b[onedem];
			gen se1_lo_demA_Pre5 = _b[onedem]-1.96*_se[onedem];
			gen se1_hi_demA_Pre5 = _b[onedem]+1.96*_se[onedem];
			gen b1_2dem_Pre5  = _b[jointdem];
			gen se1_lo_2dem_Pre5 = _b[jointdem]-1.96*_se[jointdem];
			gen se1_hi_2dem_Pre5 = _b[jointdem]+1.96*_se[jointdem];
			/* Generate placement for graph */
			gen Chg1_pre5=1;
			gen Chg1_pre4=2;
			gen Chg1_pre3=3;
			gen Chg1_pre2=4;
			gen Chg1_pre1=5;
			gen Chg1_post0=6;
			gen Chg1_post1=7;
			gen Chg1_post2=8;
			gen Chg1_post3=9;
			gen Chg1_post4=10;
			gen Chg1_post5=11;
			
			/* Collapse key variables */
			collapse se* b* Chg*;
			
			/* Graph each period for t-3 and DemA*/
			twoway rcap se3_lo_demA_Pre5 se3_hi_demA_Pre5 Chg3_pre5, lcolor(black) 
				|| rcap se3_lo_demA_Pre4 se3_hi_demA_Pre4 Chg3_pre4, lcolor(black)
				|| rcap se3_lo_demA_Pre3 se3_hi_demA_Pre3 Chg3_pre3, lcolor(black)
				|| rcap se3_lo_demA_Pre2 se3_hi_demA_Pre2 Chg3_pre2, lcolor(black)
				|| rcap se3_lo_demA_Pre1 se3_hi_demA_Pre1 Chg3_pre1, lcolor(black)
				|| rcap se3_lo_demA_Post0 se3_hi_demA_Post0 Chg3_post0, lcolor(black)
				|| rcap se3_lo_demA_Post1 se3_hi_demA_Post1 Chg3_post1, lcolor(black)
				|| rcap se3_lo_demA_Post2 se3_hi_demA_Post2 Chg3_post2, lcolor(black)
				|| rcap se3_lo_demA_Post3 se3_hi_demA_Post3 Chg3_post3, lcolor(black)
				|| rcap se3_lo_demA_Post4 se3_hi_demA_Post4 Chg3_post4, lcolor(black)
				|| rcap se3_lo_demA_Post5 se3_hi_demA_Post5 Chg3_post5, lcolor(black)
				|| scatter b3_demA_Pre5 Chg3_pre5, msymbol(circle) mcolor(black) 
				|| scatter b3_demA_Pre4 Chg3_pre4, msymbol(circle) mcolor(black) 
				|| scatter b3_demA_Pre3 Chg3_pre3, msymbol(circle) mcolor(black)
				|| scatter b3_demA_Pre2 Chg3_pre2, msymbol(circle) mcolor(black) 	
				|| scatter b3_demA_Pre1 Chg3_pre1, msymbol(circle) mcolor(black) 
				|| scatter b3_demA_Post0 Chg3_post0, msymbol(circle) mcolor(black) 
				|| scatter b3_demA_Post1 Chg3_post1, msymbol(circle) mcolor(black) 
				|| scatter b3_demA_Post2 Chg3_post2, msymbol(circle) mcolor(black) 
				|| scatter b3_demA_Post3 Chg3_post3, msymbol(circle) mcolor(black) 
				|| scatter b3_demA_Post4 Chg3_post4, msymbol(circle) mcolor(black) 
				|| scatter b3_demA_Post5 Chg3_post5, msymbol(circle) mcolor(black) 
				ytitle("") yline(0) legend(off) scheme(s1mono) t2("Democracy A") xline(5.9, lpattern(shortdash))
				xlabel(0.5 " " 1 "T-5" 2 "T-4" 3 "T-3" 4 "T-2" 5 "T-1" 6 "T" 7 "T+1" 8 "T+2" 9 "T+3" 10 "T+4" 11 "T+5" 11.5 " ", noticks)
				xticks(1(1)11);
					graph save threat3_demA, replace;
					
			/* Graph each period for t-3 and jointDem*/
			twoway rcap se3_lo_2dem_Pre5 se3_hi_2dem_Pre5 Chg3_pre5, lcolor(black) 
				|| rcap se3_lo_2dem_Pre4 se3_hi_2dem_Pre4 Chg3_pre4, lcolor(black)
				|| rcap se3_lo_2dem_Pre3 se3_hi_2dem_Pre3 Chg3_pre3, lcolor(black)
				|| rcap se3_lo_2dem_Pre2 se3_hi_2dem_Pre2 Chg3_pre2, lcolor(black)
				|| rcap se3_lo_2dem_Pre1 se3_hi_2dem_Pre1 Chg3_pre1, lcolor(black)
				|| rcap se3_lo_2dem_Post0 se3_hi_2dem_Post0 Chg3_post0, lcolor(black)
				|| rcap se3_lo_2dem_Post1 se3_hi_2dem_Post1 Chg3_post1, lcolor(black)
				|| rcap se3_lo_2dem_Post2 se3_hi_2dem_Post2 Chg3_post2, lcolor(black)
				|| rcap se3_lo_2dem_Post3 se3_hi_2dem_Post3 Chg3_post3, lcolor(black)
				|| rcap se3_lo_2dem_Post4 se3_hi_2dem_Post4 Chg3_post4, lcolor(black)
				|| rcap se3_lo_2dem_Post5 se3_hi_2dem_Post5 Chg3_post5, lcolor(black)
				|| scatter b3_2dem_Pre5 Chg3_pre5, msymbol(circle) mcolor(black) 
				|| scatter b3_2dem_Pre4 Chg3_pre4, msymbol(circle) mcolor(black) 
				|| scatter b3_2dem_Pre3 Chg3_pre3, msymbol(circle) mcolor(black)
				|| scatter b3_2dem_Pre2 Chg3_pre2, msymbol(circle) mcolor(black) 	
				|| scatter b3_2dem_Pre1 Chg3_pre1, msymbol(circle) mcolor(black) 
				|| scatter b3_2dem_Post0 Chg3_post0, msymbol(circle) mcolor(black) 
				|| scatter b3_2dem_Post1 Chg3_post1, msymbol(circle) mcolor(black) 
				|| scatter b3_2dem_Post2 Chg3_post2, msymbol(circle) mcolor(black) 
				|| scatter b3_2dem_Post3 Chg3_post3, msymbol(circle) mcolor(black) 
				|| scatter b3_2dem_Post4 Chg3_post4, msymbol(circle) mcolor(black) 
				|| scatter b3_2dem_Post5 Chg3_post5, msymbol(circle) mcolor(black) 
				ytitle("") yline(0) legend(off) scheme(s1mono) t2("Joint Democracy") xline(5.9, lpattern(shortdash))
				xlabel(0.5 " " 1 "T-5" 2 "T-4" 3 "T-3" 4 "T-2" 5 "T-1" 6 "T" 7 "T+1" 8 "T+2" 9 "T+3" 10 "T+4" 11 "T+5" 11.5 " ", noticks)
				xticks(1(1)11);
					graph save threat3_2dem, replace;

			graph combine threat3_demA.gph threat3_2dem.gph, col(2)   scheme(s1mono) title("{&Delta} Threat, t-3");
				graph save threat3, replace;

			/* Graph each period for t-2 and DemA*/
			twoway rcap se2_lo_demA_Pre5 se2_hi_demA_Pre5 Chg2_pre5, lcolor(black) 
				|| rcap se2_lo_demA_Pre4 se2_hi_demA_Pre4 Chg2_pre4, lcolor(black)
				|| rcap se2_lo_demA_Pre3 se2_hi_demA_Pre3 Chg2_pre3, lcolor(black)
				|| rcap se2_lo_demA_Pre2 se2_hi_demA_Pre2 Chg2_pre2, lcolor(black)
				|| rcap se2_lo_demA_Pre1 se2_hi_demA_Pre1 Chg2_pre1, lcolor(black)
				|| rcap se2_lo_demA_Post0 se2_hi_demA_Post0 Chg2_post0, lcolor(black)
				|| rcap se2_lo_demA_Post1 se2_hi_demA_Post1 Chg2_post1, lcolor(black)
				|| rcap se2_lo_demA_Post2 se2_hi_demA_Post2 Chg2_post2, lcolor(black)
				|| rcap se2_lo_demA_Post3 se2_hi_demA_Post3 Chg2_post3, lcolor(black)
				|| rcap se2_lo_demA_Post4 se2_hi_demA_Post4 Chg2_post4, lcolor(black)
				|| rcap se2_lo_demA_Post5 se2_hi_demA_Post5 Chg2_post5, lcolor(black)
				|| scatter b2_demA_Pre5 Chg2_pre5, msymbol(circle) mcolor(black) 
				|| scatter b2_demA_Pre4 Chg2_pre4, msymbol(circle) mcolor(black) 
				|| scatter b2_demA_Pre3 Chg2_pre3, msymbol(circle) mcolor(black)
				|| scatter b2_demA_Pre2 Chg2_pre2, msymbol(circle) mcolor(black) 	
				|| scatter b2_demA_Pre1 Chg2_pre1, msymbol(circle) mcolor(black) 
				|| scatter b2_demA_Post0 Chg2_post0, msymbol(circle) mcolor(black) 
				|| scatter b2_demA_Post1 Chg2_post1, msymbol(circle) mcolor(black) 
				|| scatter b2_demA_Post2 Chg2_post2, msymbol(circle) mcolor(black) 
				|| scatter b2_demA_Post3 Chg2_post3, msymbol(circle) mcolor(black) 
				|| scatter b2_demA_Post4 Chg2_post4, msymbol(circle) mcolor(black) 
				|| scatter b2_demA_Post5 Chg2_post5, msymbol(circle) mcolor(black) 
				ytitle("") yline(0) legend(off) scheme(s1mono) t2("Democracy A") xline(5.9, lpattern(shortdash))
				xlabel(0.5 " " 1 "T-5" 2 "T-4" 3 "T-3" 4 "T-2" 5 "T-1" 6 "T" 7 "T+1" 8 "T+2" 9 "T+3" 10 "T+4" 11 "T+5" 11.5 " ", noticks)
				xticks(1(1)11);
					graph save threat2_demA, replace;
					
			/* Graph each period for t-2 and jointDem*/
			twoway rcap se2_lo_2dem_Pre5 se2_hi_2dem_Pre5 Chg2_pre5, lcolor(black) 
				|| rcap se2_lo_2dem_Pre4 se2_hi_2dem_Pre4 Chg2_pre4, lcolor(black)
				|| rcap se2_lo_2dem_Pre3 se2_hi_2dem_Pre3 Chg2_pre3, lcolor(black)
				|| rcap se2_lo_2dem_Pre2 se2_hi_2dem_Pre2 Chg2_pre2, lcolor(black)
				|| rcap se2_lo_2dem_Pre1 se2_hi_2dem_Pre1 Chg2_pre1, lcolor(black)
				|| rcap se2_lo_2dem_Post0 se2_hi_2dem_Post0 Chg2_post0, lcolor(black)
				|| rcap se2_lo_2dem_Post1 se2_hi_2dem_Post1 Chg2_post1, lcolor(black)
				|| rcap se2_lo_2dem_Post2 se2_hi_2dem_Post2 Chg2_post2, lcolor(black)
				|| rcap se2_lo_2dem_Post3 se2_hi_2dem_Post3 Chg2_post3, lcolor(black)
				|| rcap se2_lo_2dem_Post4 se2_hi_2dem_Post4 Chg2_post4, lcolor(black)
				|| rcap se2_lo_2dem_Post5 se2_hi_2dem_Post5 Chg2_post5, lcolor(black)
				|| scatter b2_2dem_Pre5 Chg2_pre5, msymbol(circle) mcolor(black) 
				|| scatter b2_2dem_Pre4 Chg2_pre4, msymbol(circle) mcolor(black) 
				|| scatter b2_2dem_Pre3 Chg2_pre3, msymbol(circle) mcolor(black)
				|| scatter b2_2dem_Pre2 Chg2_pre2, msymbol(circle) mcolor(black) 	
				|| scatter b2_2dem_Pre1 Chg2_pre1, msymbol(circle) mcolor(black) 
				|| scatter b2_2dem_Post0 Chg2_post0, msymbol(circle) mcolor(black) 
				|| scatter b2_2dem_Post1 Chg2_post1, msymbol(circle) mcolor(black) 
				|| scatter b2_2dem_Post2 Chg2_post2, msymbol(circle) mcolor(black) 
				|| scatter b2_2dem_Post3 Chg2_post3, msymbol(circle) mcolor(black) 
				|| scatter b2_2dem_Post4 Chg2_post4, msymbol(circle) mcolor(black) 
				|| scatter b2_2dem_Post5 Chg2_post5, msymbol(circle) mcolor(black) 
				ytitle("") yline(0) legend(off) scheme(s1mono) t2("Joint Democracy") xline(5.9, lpattern(shortdash))
				xlabel(0.5 " " 1 "T-5" 2 "T-4" 3 "T-3" 4 "T-2" 5 "T-1" 6 "T" 7 "T+1" 8 "T+2" 9 "T+3" 10 "T+4" 11 "T+5" 11.5 " ", noticks)
				xticks(1(1)11);
					graph save threat2_2dem, replace;

			graph combine threat2_demA.gph threat2_2dem.gph, col(2)   scheme(s1mono) title("{&Delta} Threat, t-2");
				graph save threat2, replace;

			/* Graph each period for t-1 and DemA*/
			twoway rcap se1_lo_demA_Pre5 se1_hi_demA_Pre5 Chg1_pre5, lcolor(black) 
				|| rcap se1_lo_demA_Pre4 se1_hi_demA_Pre4 Chg1_pre4, lcolor(black)
				|| rcap se1_lo_demA_Pre3 se1_hi_demA_Pre3 Chg1_pre3, lcolor(black)
				|| rcap se1_lo_demA_Pre2 se1_hi_demA_Pre2 Chg1_pre2, lcolor(black)
				|| rcap se1_lo_demA_Pre1 se1_hi_demA_Pre1 Chg1_pre1, lcolor(black)
				|| rcap se1_lo_demA_Post0 se1_hi_demA_Post0 Chg1_post0, lcolor(black)
				|| rcap se1_lo_demA_Post1 se1_hi_demA_Post1 Chg1_post1, lcolor(black)
				|| rcap se1_lo_demA_Post2 se1_hi_demA_Post2 Chg1_post2, lcolor(black)
				|| rcap se1_lo_demA_Post3 se1_hi_demA_Post3 Chg1_post3, lcolor(black)
				|| rcap se1_lo_demA_Post4 se1_hi_demA_Post4 Chg1_post4, lcolor(black)
				|| rcap se1_lo_demA_Post5 se1_hi_demA_Post5 Chg1_post5, lcolor(black)
				|| scatter b1_demA_Pre5 Chg1_pre5, msymbol(circle) mcolor(black) 
				|| scatter b1_demA_Pre4 Chg1_pre4, msymbol(circle) mcolor(black) 
				|| scatter b1_demA_Pre3 Chg1_pre3, msymbol(circle) mcolor(black)
				|| scatter b1_demA_Pre2 Chg1_pre2, msymbol(circle) mcolor(black) 	
				|| scatter b1_demA_Pre1 Chg1_pre1, msymbol(circle) mcolor(black) 
				|| scatter b1_demA_Post0 Chg1_post0, msymbol(circle) mcolor(black) 
				|| scatter b1_demA_Post1 Chg1_post1, msymbol(circle) mcolor(black) 
				|| scatter b1_demA_Post2 Chg1_post2, msymbol(circle) mcolor(black) 
				|| scatter b1_demA_Post3 Chg1_post3, msymbol(circle) mcolor(black) 
				|| scatter b1_demA_Post4 Chg1_post4, msymbol(circle) mcolor(black) 
				|| scatter b1_demA_Post5 Chg1_post5, msymbol(circle) mcolor(black) 
				ytitle("") yline(0) legend(off) scheme(s1mono) t2("Democracy A") xline(5.9, lpattern(shortdash))
				xlabel(0.5 " " 1 "T-5" 2 "T-4" 3 "T-3" 4 "T-2" 5 "T-1" 6 "T" 7 "T+1" 8 "T+2" 9 "T+3" 10 "T+4" 11 "T+5" 11.5 " ", noticks)
				xticks(1(1)11);
					graph save threat1_demA, replace;
					
			/* Graph each period for t-1 and jointDem*/
			twoway rcap se1_lo_2dem_Pre5 se1_hi_2dem_Pre5 Chg1_pre5, lcolor(black) 
				|| rcap se1_lo_2dem_Pre4 se1_hi_2dem_Pre4 Chg1_pre4, lcolor(black)
				|| rcap se1_lo_2dem_Pre3 se1_hi_2dem_Pre3 Chg1_pre3, lcolor(black)
				|| rcap se1_lo_2dem_Pre2 se1_hi_2dem_Pre2 Chg1_pre2, lcolor(black)
				|| rcap se1_lo_2dem_Pre1 se1_hi_2dem_Pre1 Chg1_pre1, lcolor(black)
				|| rcap se1_lo_2dem_Post0 se1_hi_2dem_Post0 Chg1_post0, lcolor(black)
				|| rcap se1_lo_2dem_Post1 se1_hi_2dem_Post1 Chg1_post1, lcolor(black)
				|| rcap se1_lo_2dem_Post2 se1_hi_2dem_Post2 Chg1_post2, lcolor(black)
				|| rcap se1_lo_2dem_Post3 se1_hi_2dem_Post3 Chg1_post3, lcolor(black)
				|| rcap se1_lo_2dem_Post4 se1_hi_2dem_Post4 Chg1_post4, lcolor(black)
				|| rcap se1_lo_2dem_Post5 se1_hi_2dem_Post5 Chg1_post5, lcolor(black)
				|| scatter b1_2dem_Pre5 Chg1_pre5, msymbol(circle) mcolor(black) 
				|| scatter b1_2dem_Pre4 Chg1_pre4, msymbol(circle) mcolor(black) 
				|| scatter b1_2dem_Pre3 Chg1_pre3, msymbol(circle) mcolor(black)
				|| scatter b1_2dem_Pre2 Chg1_pre2, msymbol(circle) mcolor(black) 	
				|| scatter b1_2dem_Pre1 Chg1_pre1, msymbol(circle) mcolor(black) 
				|| scatter b1_2dem_Post0 Chg1_post0, msymbol(circle) mcolor(black) 
				|| scatter b1_2dem_Post1 Chg1_post1, msymbol(circle) mcolor(black) 
				|| scatter b1_2dem_Post2 Chg1_post2, msymbol(circle) mcolor(black) 
				|| scatter b1_2dem_Post3 Chg1_post3, msymbol(circle) mcolor(black) 
				|| scatter b1_2dem_Post4 Chg1_post4, msymbol(circle) mcolor(black) 
				|| scatter b1_2dem_Post5 Chg1_post5, msymbol(circle) mcolor(black) 
				ytitle("") yline(0) legend(off) scheme(s1mono) t2("Joint Democracy") xline(5.9, lpattern(shortdash))
				xlabel(0.5 " " 1 "T-5" 2 "T-4" 3 "T-3" 4 "T-2" 5 "T-1" 6 "T" 7 "T+1" 8 "T+2" 9 "T+3" 10 "T+4" 11 "T+5" 11.5 " ", noticks)
				xticks(1(1)11);
					graph save threat1_2dem, replace;

				graph combine threat1_demA.gph threat1_2dem.gph, col(2)   scheme(s1mono) title("{&Delta} Threat, t-1");
					graph save threat1, replace;

	/* Figure A.2 */
	graph combine threat3.gph threat2.gph threat1.gph, col(1) scheme(s1mono);
	graph save figA2, replace;

/* Table A.11 */
		/* Democracy */
		/* Load data set*/
		use NiemanGibler_demdiff_ally, clear;
		/* Replicate Leeds et al to identify shared observations*/
		logit violate demdA wcchangedA capcheither20 regchangeeither chTEAltorBmt alformA  
			asymmetry  nomicoop treaty milinst time timesquare timecube
			, robust cluster(atopidphase);
		gen LMV=1 if e(sample)==1; 
		/* Run Relevancy equation to identify shared observations*/
		logit violate Threat_combine num_dems_on_border  num_borders  cold 
			Lenergy_pop Oil_exp  maj_pow kgd_riv , cluster(atopidphase) robust;
		gen GT=1 if e(sample)==1;
		matrix B1 = e(b);
		/*Idenity cases in both samples*/
		gen GL=0;
		replace GL=1 if GT==1 & LMV==1;
		/* Rerun on Leeds et al on same sample (GL) to get starting values */
		logit violate demdA wcchangedA capcheither20 regchangeeither chTEAltorBmt alformA  
			asymmetry  nomicoop treaty milinst time timesquare timecube if GL==1
			, robust cluster (atopidphase);
		matrix B2 = e(b);
		/* Estimate Initial Model Analysis with Threat as mediator and Dem as Treatment*/
		ml model lf spl_lf  
			(B: violate = demdA wcchangedA capcheither20 regchangeeither chTEAltorBmt alformA  
			asymmetry  nomicoop treaty milinst time timesquare timecube)
			(G: violate = Threat_combine num_dems_on_border   num_borders  
			 cold Lenergy_pop Oil_exp maj_pow kgd_riv  ) if GL==1 
			, robust cluster(atopidphase) tech(nr);
		ml init B2 B1 , copy ;
		ml search;
		ml maximize, diff;
			estout , cells("b(star fmt(3)) se(par fmt(3))" ) starlevels(* .05 ) style(tex);
			estimates store DemA;
			matrix Total_B = e(b);
			matrix Total_V = e(V);
		/* Mediator equation */	
		reg Threat_combine demdA wcchangedA capcheither20 regchangeeither chTEAltorBmt alformA  
			asymmetry  nomicoop treaty milinst time timesquare timecube
			 num_dems_on_border   num_borders  
			 cold Lenergy_pop Oil_exp maj_pow kgd_riv if GL==1, robust cluster (atopidphase);
		estimates store Mediator;
		matrix Mediator_B = e(b);
		matrix Mediator_V = e(V);
		
		/*Quantities of Interest using Algorithm 1 from Imai, Keele, and Tingley 2010, p 317. */
		/* Step 1: Fit models, done above */
		/* Step 2: Simulate model parameters from sampling distribution */
		keep if GL==1; 
		expand 1000;
		drawnorm m1 m2 m3 m4 m5 m6 m7 m8 m9 m10 m11 m12 m13 m14 m15 m16 m17 m18 m19 m20 m0, mean(Mediator_B) cov(Mediator_V);
		drawnorm p1 p2 p3 p4 p5 p6 p7 p8 p9 p10 p11 p12 p13 p01 p14 p15 p16 p17 p18 p19 p20 p21 p02, mean(Total_B) cov(Total_V);
		/* Step 3 Part 1: Simulate Potential values of Mediator*/
		gen Sim_Mediator_T= m0 + m1*demdA + m2*wcchangedA + m3*capcheither20 + m4*regchangeeither + m5*chTEAltorBmt 
			+ m6*alformA + m7*asymmetry + m8*nomicoop + m9*treaty + m10*milinst + m11*time + m12*timesquare + m13*timecube
			+ m14*num_dems_on_border + m15*num_borders + m16*cold + m17*Lenergy_pop + m18*Oil_exp +m19*maj_pow + m20*kgd_riv;
		gen Sim_Mediator_C= m0 +  m2*wcchangedA + m3*capcheither20 + m4*regchangeeither + m5*chTEAltorBmt 
			+ m6*alformA + m7*asymmetry + m8*nomicoop + m9*treaty + m10*milinst + m11*time + m12*timesquare + m13*timecube
			+ m14*num_dems_on_border + m15*num_borders + m16*cold + m17*Lenergy_pop + m18*Oil_exp +m19*maj_pow + m20*kgd_riv;
		/* Step 3 Part 2: Simulate Potential Outcomes with simulated mediator values */
		gen Sim_PO_T_MT = (1/(1+exp(-(p1*demdA + p2*wcchangedA + p3*capcheither20 + p4*regchangeeither + p5*chTEAltorBmt + p6*alformA  
			+ p7*asymmetry + p8*nomicoop + p9*treaty + p10*milinst + p11*time + p12*timesquare + p13*timecube + p01))))*
			(1/(1+exp(-(p14*Sim_Mediator_T + p15*num_dems_on_border +p16*num_borders + p17*cold + p18*Lenergy_pop 
			+ p19*Oil_exp + p20*maj_pow +p21*kgd_riv+ p02)))); 
		gen Sim_PO_T_MC = (1/(1+exp(-(p1*demdA + p2*wcchangedA + p3*capcheither20 + p4*regchangeeither + p5*chTEAltorBmt + p6*alformA  
			+ p7*asymmetry + p8*nomicoop + p9*treaty + p10*milinst + p11*time + p12*timesquare + p13*timecube + p01))))*
			(1/(1+exp(-(p14*Sim_Mediator_C + p15*num_dems_on_border +p16*num_borders + p17*cold + p18*Lenergy_pop 
			+ p19*Oil_exp + p20*maj_pow +p21*kgd_riv+ p02)))); 
		gen Sim_PO_C_MT = (1/(1+exp(-( p2*wcchangedA + p3*capcheither20 + p4*regchangeeither + p5*chTEAltorBmt + p6*alformA  
			+ p7*asymmetry + p8*nomicoop + p9*treaty + p10*milinst + p11*time + p12*timesquare + p13*timecube + p01))))*
			(1/(1+exp(-(p14*Sim_Mediator_T + p15*num_dems_on_border +p16*num_borders + p17*cold + p18*Lenergy_pop 
			+ p19*Oil_exp + p20*maj_pow +p21*kgd_riv+ p02)))); 
		gen Sim_PO_C_MC = (1/(1+exp(-( p2*wcchangedA + p3*capcheither20 + p4*regchangeeither + p5*chTEAltorBmt + p6*alformA  
			+ p7*asymmetry + p8*nomicoop + p9*treaty + p10*milinst + p11*time + p12*timesquare + p13*timecube + p01))))*
			(1/(1+exp(-(p14*Sim_Mediator_C + p15*num_dems_on_border +p16*num_borders + p17*cold + p18*Lenergy_pop 
			+ p19*Oil_exp + p20*maj_pow +p21*kgd_riv+ p02)))); 
		/* Step 4: Calculate quantities of interest--Mediation, Direct, Total Effects for T,C conditions */			
		gen C_Med_T = Sim_PO_T_MT - Sim_PO_T_MC;
		gen C_Med_C = Sim_PO_C_MT - Sim_PO_C_MC;
		gen C_Med_Avg = (C_Med_T+C_Med_C)/2;
		gen C_Dir_T = Sim_PO_T_MT - Sim_PO_C_MT;
		gen C_Dir_C = Sim_PO_T_MC - Sim_PO_C_MC;
		gen C_Dir_Avg = (C_Dir_T+C_Dir_C)/2;
		gen C_Tot_T = C_Med_T + C_Dir_T;
		gen C_Tot_C = C_Med_C + C_Dir_C;
		gen C_Tot_Avg = (C_Tot_T+C_Tot_C)/2;
		collapse (mean) C_Med_T C_Med_C C_Med_Avg C_Dir_T C_Dir_C C_Dir_Avg C_Tot_T C_Tot_C C_Tot_Avg
				 (sd) C_Med_T_sd=C_Med_T C_Med_C_sd=C_Med_C C_Med_Avg_sd=C_Med_Avg
					  C_Dir_T_sd=C_Dir_T C_Dir_C_sd=C_Dir_C C_Dir_Avg_sd=C_Dir_Avg
					  C_Tot_T_sd=C_Tot_T C_Tot_C_sd=C_Tot_C C_Tot_Avg_sd=C_Tot_Avg;
		sum C_Med_T C_Med_T_sd;
		sum C_Med_C C_Med_C_sd;
		sum C_Med_Avg C_Med_Avg_sd;
		sum C_Dir_T C_Dir_T_sd;
		sum C_Dir_C C_Dir_C_sd;
		sum C_Dir_Avg C_Dir_Avg_sd;		
		sum C_Tot_T C_Tot_T_sd;
		sum C_Tot_C C_Tot_C_sd;
		sum C_Tot_Avg C_Tot_Avg_sd;			
		
		/* Democratic Allies*/	
		/* Load data set*/
		use NiemanGibler_demdiff_ally, clear;
		/* Generate indicator whether an alliance consists of two democratic members */
		gen DemAlliance = 0;
			replace DemAlliance=1 if demdA==1 & demdB==1;
		/* Replicate Leeds et al to identify shared observations*/
		logit violate wcchangedA capcheither20 regchangeeither chTEAltorBmt alformA demdA 
			asymmetry  nomicoop treaty milinst time timesquare timecube
			, robust cluster(atopidphase);
		gen LMV=1 if e(sample)==1; 
		/* Run Relevancy equation to identify shared observations*/
		logit violate Threat_combine DemAlliance num_dems_on_border  num_borders  cold 
			Lenergy_pop Oil_exp  maj_pow kgd_riv , cluster(atopidphase) robust;
		gen GT=1 if e(sample)==1;
		matrix B1 = e(b);
		/*Idenity cases in both samples*/
		gen GL=0;
		replace GL=1 if GT==1 & LMV==1;
		/* Rerun on Leeds et al on same sample (GL) to get starting values */
		logit violate wcchangedA capcheither20 regchangeeither chTEAltorBmt alformA demdA 
			asymmetry  nomicoop treaty milinst time timesquare timecube if GL==1
			, robust cluster (atopidphase);
		matrix B2 = e(b);
		/* Estimate Initial Model Analysis with Threat as mediator and Dem Ally as Treatment*/
		ml model lf spl_lf  
			(B: violate = wcchangedA capcheither20 regchangeeither chTEAltorBmt alformA demdA 
			asymmetry  nomicoop treaty milinst time timesquare timecube)
			(G: violate = Threat_combine DemAlliance num_dems_on_border   num_borders  
			 cold Lenergy_pop Oil_exp maj_pow kgd_riv  ) if GL==1 
			, robust cluster(atopidphase) tech(nr);
		ml init B2 B1 , copy ;
		ml search;
		ml maximize, diff;	
		estout , cells("b(star fmt(3)) se(par fmt(3))" ) starlevels(* .05 ) style(tex);
			estimates store DemAllies;
			matrix Total_B = e(b);
			matrix Total_V = e(V);
		/* Mediator equation */	
		reg Threat_combine DemAlliance wcchangedA capcheither20 regchangeeither chTEAltorBmt alformA demdA 
			asymmetry  nomicoop treaty milinst time timesquare timecube
			 num_dems_on_border   num_borders  
			 cold Lenergy_pop Oil_exp maj_pow kgd_riv if GL==1, robust cluster (atopidphase);
		estimates store Mediator;
		matrix Mediator_B = e(b);
		matrix Mediator_V = e(V);
		/*Quantities of Interest using Algorithm 1 from Imai, Keele, and Tingley 2010, p 317. */
		/* Step 1: Fit models, done above */
		/* Step 2: Simulate model parameters from sampling distribution */
		keep if GL==1; 
			expand 1000;
			drawnorm m1 m2 m3 m4 m5 m6 m7 m8 m9 m10 m11 m12 m13 m14 m15 m16 m17 m18 m19 m20 m21 m0, mean(Mediator_B) cov(Mediator_V);
			drawnorm p1 p2 p3 p4 p5 p6 p7 p8 p9 p10 p11 p12 p13 p01 p14 p15 p16 p17 p18 p19 p20 p21 p22 p02, mean(Total_B) cov(Total_V);
			/* Step 3 Part 1: Simulate Potential values of Mediator*/
			gen Sim_Mediator_T= m0 + m1*DemAlliance + m2*wcchangedA + m3*capcheither20 + m4*regchangeeither + m5*chTEAltorBmt 
				+ m6*alformA + m7*demdA + m8*asymmetry + m9*nomicoop + m10*treaty + m11*milinst + m12*time + m13*timesquare + m14*timecube
				+ m15*num_dems_on_border + m16*num_borders + m17*cold + m18*Lenergy_pop + m19*Oil_exp +m20*maj_pow + m21*kgd_riv;
			gen Sim_Mediator_C= m0 +  m2*wcchangedA + m3*capcheither20 + m4*regchangeeither + m5*chTEAltorBmt 
				+ m6*alformA + m7*demdA + m8*asymmetry + m9*nomicoop + m10*treaty + m11*milinst + m12*time + m13*timesquare + m14*timecube
				+ m15*num_dems_on_border + m16*num_borders + m17*cold + m18*Lenergy_pop + m19*Oil_exp +m20*maj_pow + m21*kgd_riv;
		/* Step 3 Part 2: Simulate Potential Outcomes with simulated mediator values */
			gen Sim_PO_T_MT = (1/(1+exp(-(p1*wcchangedA + p2*capcheither20 + p3*regchangeeither + p4*chTEAltorBmt + p5*alformA + p6*demdA 
				+ p7*asymmetry + p8*nomicoop + p9*treaty + p10*milinst + p11*time + p12*timesquare + p13*timecube + p01))))*
				(1/(1+exp(-(p14*Sim_Mediator_T + p15*DemAlliance + p16*num_dems_on_border +p17*num_borders + p18*cold + p19*Lenergy_pop 
				+ p20*Oil_exp + p21*maj_pow +p22*kgd_riv+ p02)))); 
			gen Sim_PO_T_MC = (1/(1+exp(-(p1*wcchangedA + p2*capcheither20 + p3*regchangeeither + p4*chTEAltorBmt + p5*alformA + p6*demdA 
				+ p7*asymmetry + p8*nomicoop + p9*treaty + p10*milinst + p11*time + p12*timesquare + p13*timecube + p01))))*
				(1/(1+exp(-(p14*Sim_Mediator_C + p15*DemAlliance + p16*num_dems_on_border +p17*num_borders + p18*cold + p19*Lenergy_pop 
				+ p20*Oil_exp + p21*maj_pow +p22*kgd_riv+ p02)))); 
			gen Sim_PO_C_MT = (1/(1+exp(-( p2*capcheither20 + p3*regchangeeither + p4*chTEAltorBmt + p5*alformA + p6*demdA 
				+ p7*asymmetry + p8*nomicoop + p9*treaty + p10*milinst + p11*time + p12*timesquare + p13*timecube + p01))))*
				(1/(1+exp(-(p14*Sim_Mediator_T + p15*DemAlliance + p16*num_dems_on_border +p17*num_borders + p18*cold + p19*Lenergy_pop 
				+ p20*Oil_exp + p21*maj_pow +p22*kgd_riv+ p02)))); 
			gen Sim_PO_C_MC = (1/(1+exp(-( p2*capcheither20 + p3*regchangeeither + p4*chTEAltorBmt + p5*alformA + p6*demdA 
				+ p7*asymmetry + p8*nomicoop + p9*treaty + p10*milinst + p11*time + p12*timesquare + p13*timecube + p01))))*
				(1/(1+exp(-(p14*Sim_Mediator_C + p15*DemAlliance + p16*num_dems_on_border +p17*num_borders + p18*cold + p19*Lenergy_pop 
				+ p20*Oil_exp + p21*maj_pow +p22*kgd_riv+ p02))));
		/* Step 4: Calculate quantities of interest--Mediation, Direct, Total Effects for T,C conditions */			
			gen C_Med_T = Sim_PO_T_MT - Sim_PO_T_MC;
			gen C_Med_C = Sim_PO_C_MT - Sim_PO_C_MC;
			gen C_Med_Avg = (C_Med_T+C_Med_C)/2;
			gen C_Dir_T = Sim_PO_T_MT - Sim_PO_C_MT;
			gen C_Dir_C = Sim_PO_T_MC - Sim_PO_C_MC;
			gen C_Dir_Avg = (C_Dir_T+C_Dir_C)/2;
			gen C_Tot_T = C_Med_T + C_Dir_T;
			gen C_Tot_C = C_Med_C + C_Dir_C;
			gen C_Tot_Avg = (C_Tot_T+C_Tot_C)/2;
			collapse (mean) C_Med_T C_Med_C C_Med_Avg C_Dir_T C_Dir_C C_Dir_Avg C_Tot_T C_Tot_C C_Tot_Avg
					 (sd) C_Med_T_sd=C_Med_T C_Med_C_sd=C_Med_C C_Med_Avg_sd=C_Med_Avg
						  C_Dir_T_sd=C_Dir_T C_Dir_C_sd=C_Dir_C C_Dir_Avg_sd=C_Dir_Avg
						  C_Tot_T_sd=C_Tot_T C_Tot_C_sd=C_Tot_C C_Tot_Avg_sd=C_Tot_Avg;
			sum C_Med_T C_Med_T_sd;
			sum C_Med_C C_Med_C_sd;
			sum C_Med_Avg C_Med_Avg_sd;
			sum C_Dir_T C_Dir_T_sd;
			sum C_Dir_C C_Dir_C_sd;
			sum C_Dir_Avg C_Dir_Avg_sd;		
			sum C_Tot_T C_Tot_T_sd;
			sum C_Tot_C C_Tot_C_sd;
			sum C_Tot_Avg C_Tot_Avg_sd;	
		
	/* Table A.11 Results 
	(quantities of interest for lower part of table from Step 4 for Dem and Dem Allies above)*/
		estout DemA DemAllies, cells("b(star fmt(3)) se(par fmt(3))" ) starlevels(* .05 ) style(tex);		
		
/* Table A.12 */
		/* Load data set and organize data */
		use NiemanGibler_demdiff_ally, clear;
		sort year ccode1 ccode2;
	/* Table A.12 Results*/
	list ccode1 ccode2 year violate Threat_combine wcchangedA demdA if violate==1;	
	
/* Table A.13 */
		/* Calculate Rubin formula for uncertainty adding R-squared */
		capture program drop rubin;
		program define rubin;
			local names: colnames(Q1);
			/* MI param estimates */
			matrix Q = (Q1 + Q2 + Q3 + Q4 + Q5 + Q6 + Q7 + Q8 + Q9 + Q10)/10;
			/* MI covariances     */
			matrix W = (W1 + W2 + W3 + W4 + W5 + W6 + W7 + W8 + W9 + W10)/10;
			/* MI R-squared */
			matrix r2 = (r21 + r22 + r23 + r24 + r25 + r26 + r27 + r28 + r29 + r210)/10;
			local  k = colsof(Q);
			forvalues i=1/10 {;
			  matrix QQ=Q`i'-Q;
			  if `i'==1 {;
				matrix B=(QQ)'*QQ;
				};
			  else matrix B = B + (QQ)'*QQ;
			};
			/* Covariance adjustment */
			matrix B=B/(10-1);		
			/* Final covariance   */
			matrix V=W+(1+1/10)*B;				
			ereturn post Q V ;
		end;
		/* Load data */
		use NiemanGibler_demdiff_trade, clear;
		/* replicate original findings */
		reg ltrade bothin onein gsp ldist lrgdp lrgdppc regional custrict comlang border landl island lareap comcol colony curcol comctry i.year, robust cluster(dyad);
		estimates store rose;		
		/* Adding Threat (main model)*/
		forvalues i=1/10 {;
		areg ltrade BothDem_Threatmin`i' OneDem_Threatmin`i' BothDem OneDem threat_min`i' bothin onein gsp ldist lrgdp lrgdppc regional custrict comlang border landl island lareap comcol colony curcol comctry , a(year) robust cluster(dyad);
			matrix Q`i' = e(b);
			matrix W`i' = e(V);
			matrix  r2`i' = e(r2);
			nois _dots `i' 0;
		};
		rubin;
			estout , cells(b(star fmt(3)) se(par fmt(3)) ) starlevels(* .05) style(tex);
			matrix B = e(b);
			matrix V = e(V);
			matrix rename  r2 r2_threat;
			estimates store threat;
	
		/*rerun using combined scores for interaction figures*/
		areg ltrade BothDem_Threatmin OneDem_Threatmin BothDem OneDem threat_min bothin onein gsp ldist lrgdp lrgdppc regional custrict comlang border landl island lareap comcol colony curcol comctry , a(year) robust cluster(dyad);
			gen GN =1 if e(sample)==1;
		/* Reduced Sample */
		areg ltrade bothin onein gsp ldist lrgdp lrgdppc regional custrict comlang border landl island lareap comcol colony curcol comctry  if GN==1, a(year) robust cluster(dyad);
		estimates store rose2;
		/* Adding Democracy */
		areg ltrade BothDem OneDem bothin onein gsp ldist lrgdp lrgdppc regional custrict comlang border landl island lareap comcol colony curcol comctry  if GN==1, a(year) robust cluster(dyad);
		estimates store rose_dem;	
		
	/* Table A.13 Results */	
	#delimit;
	estout rose rose2 rose_dem threat, stats(N r2)
	cells(b(fmt(3) star) se(par))
	modelwidth(8)
	starlevels(* 0.05) legend
	varwidth(15)
	drop( *year)
	style(tex) replace;
	
	matrix list r2_threat;	

/* Figure A.3 */
		/* Preserve data to restore later */
		preserve;
		/* Limit data to sample */
		keep if GN==1;
		/* Collapse to range of threat min to max */
		collapse (mean) threat=threat_min (min) threat_min=threat_min (max) threat_max=threat_min;
		/* Generate simulations */
		expand 100;
		replace threat = threat_min + (_n/_N)*(threat_max - threat_min);
	    expand 100;
		drawnorm b_BothDem_Threatmin b_OneDem_Threatmin b_BothDem b_OneDem b_threat_min b_bothin b_onein b_gsp b_ldist b_lrgdp b_lrgdppc b_regional b_custrict b_comlang b_border b_landl b_island b_lareap b_comcol b_colony b_curcol b_comctry b_cons, means(B) cov(V);
		/* Generate marginal effects */
		gen gamma1 = b_OneDem + threat*b_OneDem_Threatmin;
		gen gamma2 = b_BothDem + threat*b_BothDem_Threatmin;
		/* Get mean and SD for each value of z. */
		collapse (mean) gamma1 gamma2 (sd) sd_gamma1 = gamma1 sd_gamma2 = gamma2, by(threat);
		gen gamma1_lo = gamma1 - 1.96*sd_gamma1;
		gen gamma1_hi = gamma1 + 1.96*sd_gamma1;
		gen gamma2_lo = gamma2 - 1.96*sd_gamma2;
		gen gamma2_hi = gamma2 + 1.96*sd_gamma2; 
	/* Make figure A.3*/
		twoway lowess gamma1 threat, scheme(s1color)
			lcolor(blue) lwidth(thick)
			|| 	lowess gamma1_lo threat, 
			lcolor(ltblue ) lwidth(thin  ) 
			|| 	lowess gamma1_hi threat, 
			lcolor(ltblue ) lwidth(thin  ) 
			|| lowess gamma2 threat, 
			lcolor(red) lwidth(thick) lpattern(longdash )
			|| 	lowess gamma2_lo  threat, 
			lcolor(pink ) lwidth(thin ) lpattern(longdash )
			|| 	lowess  gamma2_hi threat, 
			lcolor(pink ) lwidth(thin ) lpattern(longdash )
			yline(0)
			xtitle("Territorial Threat (weak link)")
			legend(off);
	graph save figA3, replace asis;

/* Figure A.4*/	
	restore;	
	preserve;
		/* Limite data to sample and collapse to min/max range */
		keep if GN==1;
		collapse (mean) threat=threat_min  (median) ldist lrgdp lrgdppc lareap bothin onein gsp  regional custrict comlang border landl island comcol colony curcol comctry (min) threat_min=threat_min (max) threat_max=threat_min;
		/* Generate simulations */
		expand 100;	
		replace threat = threat_min + (_n/_N)*(threat_max - threat_min);
		expand 100;
		drawnorm b_BothDem_Threatmin b_OneDem_Threatmin b_BothDem b_OneDem b_threat_min b_bothin b_onein b_gsp b_ldist b_lrgdp b_lrgdppc b_regional b_custrict b_comlang b_border b_landl b_island b_lareap b_comcol b_colony b_curcol b_comctry b_cons, means(B) cov(V);
		gen pr_trade_both = exp(b_BothDem_Threatmin*threat + b_BothDem + b_threat_min*threat + b_bothin*bothin + b_onein*onein + b_gsp*gsp + b_ldist*ldist + b_lrgdp*lrgdp +
			b_lrgdppc*lrgdppc + b_regional*regional + b_custrict*custrict + b_comlang*comlang +  b_border*border + b_landl*landl + b_island*island + b_lareap*lareap + 
			b_comcol*comcol + b_colony*colony + b_curcol*curcol + b_comctry*comctry + b_cons);
		gen pr_trade_one = exp(b_OneDem_Threatmin*threat + b_OneDem + b_threat_min*threat + b_bothin*bothin + b_onein*onein + b_gsp*gsp + b_ldist*ldist + b_lrgdp*lrgdp +
			b_lrgdppc*lrgdppc + b_regional*regional + b_custrict*custrict + b_comlang*comlang +  b_border*border + b_landl*landl + b_island*island + b_lareap*lareap + 
			b_comcol*comcol + b_colony*colony + b_curcol*curcol + b_comctry*comctry + b_cons);
		gen pr_trade_none = exp(b_threat_min*threat + b_bothin*bothin + b_onein*onein + b_gsp*gsp + b_ldist*ldist + b_lrgdp*lrgdp +
			b_lrgdppc*lrgdppc + b_regional*regional + b_custrict*custrict + b_comlang*comlang +  b_border*border + b_landl*landl + b_island*island + b_lareap*lareap + 
			b_comcol*comcol + b_colony*colony + b_curcol*curcol + b_comctry*comctry + b_cons);
		/* Calculate quantities of interest */
		collapse (mean) pr_trade_both pr_trade_one pr_trade_none, by(threat);
	/* Make figure */
		twoway lowess pr_trade_both threat, scheme(s1color)
			lcolor(blue) 
			|| lowess pr_trade_one threat, 
			lcolor(red) lwidth(medthick) lpattern (longdash)
			|| 	lowess  pr_trade_none threat, 
			lcolor(black ) lwidth(medium ) lpattern(dash )
			ylabel(0(5000)15000)
			xtitle("Territorial Threat (weak link)")
			legend(lab(1  "Both Dem") lab(2 "One Dem") lab(3 "No Dem") symx(*.5) row(1));
			graph save figA4, replace asis;
restore;	
		
log close;
